01179 {
01180   int ip;                    /* parameter index */
01181   int nxyz;                  /* total number of voxels */
01182   int ixyz;                  /* voxel index */
01183 
01184 
01185   /*----- Initialize local variables -----*/
01186   nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01187 
01188 
01189   /*----- Allocate memory space for fim parameters -----*/
01190   *fim_params_vol  = (float **) malloc (sizeof(float *) * NBUCKETS);   
01191   MTEST(*fim_params_vol);
01192 
01193   for (ip = 0;  ip < NBUCKETS;  ip++)
01194     {
01195                                   (*fim_params_vol)[ip]  = (float *) malloc (sizeof(float) * nxyz);
01196                                   MTEST((*fim_params_vol)[ip]);    
01197                                   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01198                                  {
01199                                 (*fim_params_vol)[ip][ixyz]  = 0.0;
01200                                  }
01201     }
01202 
01203 }

01785 {
01786   const float EPSILON = 1.0e-10;
01787   float factor;             /* factor is new scale factor for sub-brick #ib */
01788 
01789 
01790   /*----- allocate memory for output sub-brick -----*/
01791   bar[ibrick]  = (short *) malloc (sizeof(short) * nxyz);
01792   MTEST (bar[ibrick]);
01793   factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01794                                       MRI_short, bar[ibrick]);
01795   
01796   if (factor < EPSILON)  factor = 0.0;
01797   else factor = 1.0 / factor;
01798   
01799 
01800   /*----- edit the sub-brick -----*/
01801   EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
01802   EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
01803 
01804   if (brick_type == FUNC_COR_TYPE)
01805     EDIT_BRICK_TO_FICO (new_dset, ibrick, nsam, nfit, nort);
01806   
01807   
01808   /*----- attach bar[ib] to be sub-brick #ibrick -----*/
01809   EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
01810 
01811 }

01344 {
01345   float * ts_array = NULL;    /* array of measured data for one voxel */
01346   float mask_val[1];          /* value of mask at current voxel */
01347 
01348   int q;                      /* number of parameters in the baseline model */
01349   int p;                      /* number of parameters in the baseline model 
01350                                  plus number of ideals */
01351   int m;                      /* parameter index */
01352   int n;                      /* data point index */
01353 
01354 
01355   matrix xdata;               /* independent variable matrix */
01356   matrix x_base;              /* extracted X matrix    for baseline model */
01357   matrix xtxinvxt_base;       /* matrix:  (1/(X'X))X'  for baseline model */
01358   matrix x_ideal[MAX_FILES];  /* extracted X matrices  for ideal models */
01359   matrix xtxinvxt_ideal[MAX_FILES];     
01360                               /* matrix:  (1/(X'X))X'  for ideal models */
01361   vector y;                   /* vector of measured data */       
01362 
01363   int ixyz;                   /* voxel index */
01364   int nxyz;                   /* number of voxels per dataset */
01365 
01366   int nt;                  /* number of images in input 3d+time dataset */
01367   int NFirst;              /* first image from input 3d+time dataset to use */
01368   int NLast;               /* last image from input 3d+time dataset to use */
01369   int N;                   /* number of usable data points */
01370 
01371   int num_ort_files;       /* number of ort time series files */
01372   int num_ideal_files;     /* number of ideal time series files */
01373   int polort;              /* degree of polynomial ort */
01374   int num_ortts;           /* number of ort time series */
01375   int num_idealts;         /* number of ideal time series */
01376   
01377   int i;                   /* data point index */
01378   int is;                  /* ideal index */
01379   int ilag;                /* time lag index */
01380   float stddev;            /* normalized parameter standard deviation */
01381   char * label;            /* string containing stat. summary of results */
01382         int xpos, ypos, zpos; 
01383    double tzero , tdelta , ts_mean , ts_slope;
01384    int   ii ,  iz,izold, nxy , nuse, use_fac, kk;
01385   float * x_bot = NULL;    /* minimum of stimulus time series */
01386   float * x_ave = NULL;    /* average of stimulus time series */
01387   float * x_top = NULL;    /* maximum of stimulus time series */
01388   int * good_list = NULL;  /* list of good (usable) time points */ 
01389   float ** rarray = NULL;  /* ranked arrays of ideal time series */
01390   float FimParams[NBUCKETS];  /* output delay parameters */
01391 
01392    float *  dtr  = NULL ;  /* will be array of detrending coeff */
01393    float *  fac  = NULL ;  /* array of input brick scaling factors */
01394         float * vox_vect;                       /* voxel time series */
01395         float *ref_ts; /*reference time series */
01396         float slp, delu, del,  xcor, xcorCoef,vts, vrvec, dtx, d0fac , d1fac , x0,x1;
01397         int  actv, opt, iposdbg;
01398         
01399         #ifdef ZDBG
01400                 xyzTOindex (option_data, &iposdbg, IPOSx,  IPOSy , IPOSz);
01401                 printf ("Debug for %d: %d, %d, %d\n\n", iposdbg, IPOSx,  IPOSy , IPOSz);
01402         #endif
01403 
01404   /*----- Initialize matrices and vectors -----*/
01405   matrix_initialize (&xdata);
01406   matrix_initialize (&x_base);
01407   matrix_initialize (&xtxinvxt_base);
01408   for (is =0;  is < MAX_FILES;  is++)
01409     {
01410       matrix_initialize (&x_ideal[is]);
01411       matrix_initialize (&xtxinvxt_ideal[is]);
01412     }
01413   vector_initialize (&y);
01414 
01415 
01416   /*----- Initialize local variables -----*/
01417   if (option_data->input1D_filename != NULL)
01418     {
01419       nxyz = 1;
01420       nt = fmri_length;
01421     }
01422   else
01423     {
01424       nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;       
01425       nt = DSET_NUM_TIMES (dset);
01426     }
01427 
01428   NFirst = option_data->NFirst;
01429   NLast = option_data->NLast;   
01430   N = option_data->N;
01431   p = option_data->p;
01432   q = option_data->q;
01433 
01434   polort          = option_data->polort;
01435   num_idealts     = option_data->num_idealts;
01436   num_ort_files   = option_data->num_ort_files;
01437   num_ideal_files = option_data->num_ideal_files;
01438 
01439 
01440   /*----- Allocate memory -----*/
01441   ts_array = (float *) malloc (sizeof(float) * nt);      MTEST (ts_array);
01442   x_bot = (float *)    malloc (sizeof(float) * p);       MTEST (x_bot);
01443   x_ave = (float *)    malloc (sizeof(float) * p);       MTEST (x_ave);
01444   x_top = (float *)    malloc (sizeof(float) * p);       MTEST (x_top);
01445   good_list = (int *)  malloc (sizeof(int) * N);         MTEST (good_list);
01446   rarray = (float **)  malloc (sizeof(float *) * num_idealts);  MTEST (rarray);
01447   vox_vect = (float *)    malloc (sizeof(float) * nt);       MTEST (vox_vect);
01448   
01449   /*----- Initialize the independent variable matrix -----*/
01450   N = init_indep_var_matrix (q, p, NFirst, N, polort, 
01451                              num_ort_files, num_ideal_files, 
01452                              ort_array, ort_list, ideal_array, ideal_list, 
01453                              x_bot, x_ave, x_top, good_list, &xdata);
01454   option_data->N = N;
01455 
01456 
01457 
01458  
01459   /*----- Initialization for the regression analysis -----*/
01460   init_delay (q, p, N, num_idealts, xdata, &x_base, 
01461                             &xtxinvxt_base, x_ideal, xtxinvxt_ideal, rarray);
01462 
01463 
01464   vector_create (N, &y);
01465 
01466    /*----- set up to find time at each voxel -----*/
01467 
01468    tdelta = dset->taxis->ttdel ;
01469    if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
01470    if( tdelta == 0.0 ) tdelta = 1.0 ;
01471    izold  = -666 ;
01472    nxy    = dset->daxes->nxx * dset->daxes->nyy ;
01473   
01474   
01475    /*--- get scaling factors for input sub-bricks ---*/
01476         nuse      = DSET_NUM_TIMES(dset) - option_data->NFirst ;
01477    fac = (float *) malloc( sizeof(float) * nuse ) ;   /* factors */ MTEST (fac);
01478    
01479 
01480    use_fac = 0 ;
01481    for( kk=0 ; kk < nuse ; kk++ ){
01482       fac[kk] = DSET_BRICK_FACTOR(dset,kk+option_data->NFirst) ;
01483       if( fac[kk] != 0.0 ) use_fac++ ;
01484       else                 fac[kk] = 1.0 ;
01485    }
01486    if( !use_fac ) FREEUP(fac) ;
01487 
01488    /*--- setup for detrending ---*/
01489 
01490    dtr = (float *) malloc( sizeof(float) * nuse ) ;MTEST (dtr);
01491    
01492 
01493    d0fac = 1.0 / nuse ;
01494    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
01495    for( kk=0 ; kk < nuse ; kk++ )
01496       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* linear trend, orthogonal to 1 */
01497 
01498 
01499   /*----- Loop over all voxels -----*/
01500   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01501     {
01502                 #ifdef ZDBG
01503                         if (ixyz == iposdbg)
01504                                 {
01505                                         printf ("\nAt voxel, checking for mask\n");
01506                                 }
01507                 #endif
01508 
01509                 indexTOxyz ( option_data , ixyz, &xpos , &ypos , &zpos);
01510       /*----- Apply mask? -----*/
01511       if (mask != NULL)
01512                         {
01513                           extract_ts_array (mask, ixyz, mask_val);
01514                           if (mask_val[0] == 0.0)  continue; 
01515                         }
01516                 
01517                 #ifdef ZDBG
01518                         if (ixyz == iposdbg)
01519                                 {
01520                                 printf ("Past the mask, extracting TS\n");
01521                                 }
01522                 #endif
01523 
01524 
01525      /*----- Extract Y-data for this voxel -----*/
01526         if (option_data->input1D_filename != NULL)
01527                         {
01528                           for (i = 0;  i < N;  i++)
01529                          vox_vect[i] = (float)fmri_data[good_list[i]+NFirst];
01530                         }
01531       else
01532                         {
01533                           extract_ts_array (dset, ixyz, ts_array);
01534                           for (i = 0;  i < N;  i++)
01535                          vox_vect[i] = (float)ts_array[good_list[i]+NFirst];
01536                         }
01537       
01538                 #ifdef ZDBG
01539                         if (ixyz == iposdbg)
01540                                 {
01541                                         printf ("TS extracted\n");
01542                                 }
01543                 #endif
01544 
01545 
01546         
01547         opt = 1; /* set to 0 for cleanup */
01548         
01549       /*** scale? ***/
01550 
01551                 #ifdef ZDBG
01552                 if (ixyz == iposdbg)
01553                         {
01554                                 printf("Before Scale\n");
01555                                 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1],  vox_vect[2], vox_vect[3], vox_vect[4]);
01556                                 /*getchar ();*/
01557                         }
01558                 #endif
01559                         
01560       if( use_fac )
01561          for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] *= fac[kk] ;
01562       
01563                 #ifdef ZDBG
01564                 if (ixyz == iposdbg)
01565                         {
01566                                 printf("Before Detrending\n");
01567                                 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1],  vox_vect[2], vox_vect[3], vox_vect[4]);
01568                                 /*getchar ();*/
01569                         }
01570                 #endif
01571 
01572                 /** compute mean and slope **/
01573 
01574       x0 = x1 = 0.0 ;
01575       for( kk=0 ; kk < nuse ; kk++ ){
01576          x0 += vox_vect[kk] ; x1 += vox_vect[kk] * dtr[kk] ;
01577       }
01578 
01579       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
01580 
01581       ts_mean  = x0 ;
01582       ts_slope = x1 / tdelta ;
01583 
01584       /** detrend? **/
01585 
01586       if( option_data->dtrnd )
01587          for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] -= (x0 + x1 * dtr[kk]) ;
01588       else
01589          for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] -= x0;
01590          
01591                 #ifdef ZDBG
01592                 if (ixyz == iposdbg)
01593                         {
01594                                 printf("After Detrending (or just zero-meaning)\n");
01595                                 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1],  vox_vect[2], vox_vect[3], vox_vect[4]);
01596                                 /*getchar ();*/
01597                         }
01598                 #endif
01599 
01600                 /* calculate the T0 and Tdelta */
01601       /** compute start time of this timeseries **/
01602 
01603       iz = ixyz / nxy ;    /* which slice am I in? */
01604 
01605       if( iz != izold ){          /* in a new slice? */
01606          tzero = THD_timeof( option_data->NFirst ,
01607                              dset->daxes->zzorg
01608                            + iz*dset->daxes->zzdel , dset->taxis ) ;
01609          izold = iz ;
01610 
01611          if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
01612      
01613           if (option_data->Dsamp == YUP) 
01614                         dtx = (float) (tzero / tdelta) - option_data->NFirst;
01615                 else
01616                         dtx = 0.0;
01617       }
01618         
01619                 #ifdef ZDBG
01620                 if (ixyz == iposdbg)
01621                         {
01622                                 printf("dtx = %f\n", dtx);
01623                         }
01624                 #endif
01625 
01626         option_data->errcode = hilbertdelay_V2 (vox_vect, option_data->rvec,option_data->ln,option_data->Nseg,option_data->Pover,opt,0,dtx,option_data->biasrem,&delu,&slp,&xcor,&xcorCoef,&vts,&vrvec);                                        /* cleanup time */
01627                 #ifdef ZDBG
01628                 if (ixyz == iposdbg)
01629                         {
01630                                 printf ("%d err code @ixyz = %d\n", option_data->errcode, ixyz);
01631                         }
01632                 #endif
01633         if (option_data->errcode == 0) /* If there are no errors, proceed */
01634                 { /*option_data->errcode == 0 inner loop */
01635                                         hunwrap (delu, (float)(option_data->fs), option_data->T, slp, option_data->wrp, option_data->unt, &del );
01636                                         
01637                                         actv = 1;                                               /* assume voxel is active */
01638         
01639                                         if (xcorCoef < option_data->co) actv = 0;                       /* determine if voxel is activated using xcorCoef  */
01640         
01641                                         if ((actv == 1) && (option_data->out == YUP))           /* if voxel is truly activated, write results to file without modifying return value */
01642                                                 {
01643                                                         indexTOxyz ( option_data , ixyz, &xpos , &ypos , &zpos);
01644                                                         fprintf (option_data->outwrite,"%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\n", ixyz , xpos , ypos , zpos ,  delu , del , xcor , xcorCoef , vts);
01645                                                         if (option_data->outts == YUP)
01646                                                                 {
01647                                                                         writets (option_data, vox_vect);        
01648                                                                 }
01649                                                 }
01650 
01651                 }/*option_data->errcode == 0 inner loop */
01652                         
01653         else if (option_data->errcode == ERROR_LONGDELAY)
01654                                 {                                       
01655                                         error_report ( option_data , ixyz);     
01656 
01657                                         del = 0.0;                                                              /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
01658                                 xcorCoef = 0.0;                                         /*  because the data might still be OK */
01659                                 xcor = 0.0;
01660 
01661                                 }
01662                         else if (option_data->errcode != 0)
01663                                 {
01664                                         error_report ( option_data , ixyz);     
01665 
01666                                         del = 0.0;                                                              /* Set all the variables to Null and set xcorCoef to an impossible value*/
01667                                 xcorCoef = NOWAYXCORCOEF;                                               
01668                                 xcor = 0.0;
01669                                 }       
01670         
01671         
01672                 FimParams[DELINDX] = del;
01673                 FimParams[COVINDX] = xcor;
01674                 FimParams[COFINDX] = xcorCoef;
01675                 FimParams[VARINDX] = vts;
01676         
01677                 #ifdef ZDBG
01678                 if (ixyz == iposdbg)
01679                         {
01680                                 printf("VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
01681                                 printf("%d\t%d\t%d\t%d\t", ixyz, xpos, ypos, zpos);
01682                                 printf("%f\t%f\t%f\t%f\t%f\n", delu, del, xcor, xcorCoef, vts);
01683                                 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1],  vox_vect[2], vox_vect[3], vox_vect[4]);
01684                                 printf("%f\t%f\t%f\t%f\t%f\n", dtx, delu, del, xcor, xcorCoef);
01685                                 /*getchar ();*/
01686                         }
01687                 #endif
01688         
01689       /*----- Save results for this voxel -----*/
01690       if (option_data->input1D_filename == NULL)
01691         save_voxel (ixyz, FimParams, fim_params_vol);
01692       
01693       
01694       
01695     }  /*----- Loop over voxels -----*/
01696   
01697 
01698    if (option_data->out == YUP)                                                                 /* close outfile and outlogfile*/
01699                                 {
01700                                         fclose (option_data->outlogfile);
01701                                         fclose (option_data->outwrite);
01702                                         if (option_data->outts  == YUP) fclose (option_data->outwritets);
01703                                 }
01704                                 else
01705                                 {
01706                                         if (option_data->outlogfile != NULL)    fclose (option_data->outlogfile);               /* close outlogfile */
01707                                 }
01708 
01709 
01710   /*----- Dispose of matrices and vectors -----*/
01711   vector_destroy (&y);
01712   for (is = 0;  is < MAX_FILES;  is++)
01713     {
01714       matrix_destroy (&xtxinvxt_ideal[is]);
01715       matrix_destroy (&x_ideal[is]);
01716     } 
01717   matrix_destroy (&xtxinvxt_base);
01718   matrix_destroy (&x_base);
01719   matrix_destroy (&xdata);
01720 
01721 
01722   /*----- Deallocate memory -----*/
01723   free (ts_array);     ts_array = NULL;
01724   free (x_bot);        x_bot = NULL;
01725   free (x_ave);        x_ave = NULL;
01726   free (x_top);        x_top = NULL;
01727   free (good_list);    good_list = NULL;
01728   free (vox_vect);              vox_vect = NULL;
01729   
01730   for (is = 0;  is < num_idealts;  is++)
01731     {
01732       free (rarray[is]);   rarray[is] = NULL;
01733     }
01734   free (rarray);       rarray = NULL;
01735   
01736 }

01026 {
01027   char message[THD_MAX_NAME];  /* error message */
01028   int is;                  /* ideal index */
01029   int num_ort_files;       /* number of ort time series files */
01030   int num_ideal_files;     /* number of ideal time series files */
01031   int num_idealts;         /* number of ideal time series */
01032   int nt;                  /* number of images in input 3d+time dataset */
01033   int NFirst;              /* first image from input 3d+time dataset to use */
01034   int NLast;               /* last image from input 3d+time dataset to use */
01035   int N;                   /* number of usable time points */
01036         int lncheck;
01037 
01038   /*----- Initialize local variables -----*/
01039   if (option_data->input1D_filename != NULL)
01040     nt = fmri_length;
01041   else
01042     nt = DSET_NUM_TIMES (dset_time);
01043 
01044   num_ort_files   = option_data->num_ort_files;
01045   num_ideal_files = option_data->num_ideal_files;
01046   num_idealts     = option_data->num_idealts;
01047 
01048 
01049   NFirst = option_data->NFirst;
01050 
01051   NLast = option_data->NLast;   
01052   if (NLast > nt-1)  NLast = nt-1;
01053   option_data->NLast = NLast;
01054 
01055   N = NLast - NFirst + 1;
01056   option_data->N = N;
01057 
01058 
01059   /*----- Check number of ideal time series -----*/
01060   if (num_idealts < 1)  FIM_error ("No ideal time series?");
01061 
01062 
01063   /*----- If mask is used, check for compatible dimensions -----*/
01064   if (mask_dset != NULL)
01065     {
01066       if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
01067            || (DSET_NY(dset_time) != DSET_NY(mask_dset))
01068            || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
01069         {
01070           sprintf (message, "%s and %s have incompatible dimensions",
01071                    option_data->input_filename, option_data->mask_filename);
01072           FIM_error (message);
01073         }
01074 
01075       if (DSET_NVALS(mask_dset) != 1 )
01076         FIM_error ("Must specify 1 sub-brick from mask dataset");
01077     }
01078 
01079 
01080 
01081  
01082   /*----- Check whether any of the output files already exist -----*/
01083   check_output_files (option_data, dset_time);
01084  
01085   /*----- Read in reference time series -----*/
01086    option_data->ln = option_data->NLast - option_data->NFirst + 1;
01087         option_data->rvec = (float *)    malloc (sizeof(float) * option_data->ln+1);       MTEST (option_data->rvec);
01088 
01089   /*------- Load Reference Time Series ------------------*/
01090   lncheck = float_file_size (option_data->ideal_filename[0]);
01091   if (lncheck != nt)
01092         {
01093                 printf("Error: Reference filename contains %d values.\n %d values were expected.\n", lncheck, nt);
01094                 exit (0);
01095         }
01096   
01097         Read_part_file_delay (option_data->rvec, option_data->ideal_filename[0], option_data->NFirst + 1,option_data->NLast + 1);  
01098 
01099          
01100   /* --- decide on the bucket name ----*/ 
01101         if (option_data->bucket_filename == NULL)
01102         {
01103                 option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
01104                 MTEST (option_data->bucket_filename);
01105                 sprintf (option_data->bucket_filename, "%s.DEL", DSET_PREFIX (dset_time));
01106                 /*make sure that prefix is OK*/
01107                 check_output_files (option_data, dset_time);
01108         }
01109         
01110   /* --- decide on the output name ----*/ 
01111         /* The log file is created no matter what */
01112         if (option_data->outname == NULL)
01113                 {
01114                 option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
01115                         MTEST (option_data->outnamelog);
01116                         sprintf (option_data->outnamelog, "%s.log", option_data->bucket_filename);
01117                 }
01118         if (option_data->out || option_data->outts)
01119         {
01120                 if (option_data->outname == NULL)
01121                 {
01122                         option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
01123                         MTEST (option_data->outname);
01124                         sprintf (option_data->outname, "%s", option_data->bucket_filename);
01125                 }
01126                 if (option_data->outts)
01127                 {
01128                         option_data->outnamets = malloc (sizeof(char)*(THD_MAX_NAME+3));
01129                         MTEST (option_data->outnamets);
01130                         sprintf (option_data->outnamets, "%s.ts", option_data->outname);
01131                 }
01132         }
01133         
01134  /* ------- Open files for writing -------------*/
01135         
01136         option_data->outlogfile = fopen (option_data->outnamelog,"w"); /* open log file regardless */
01137         
01138         if (option_data->out == YUP)                                                                    /* open outfile */
01139                                 {                                       
01140                                         option_data->outwrite = fopen (option_data->outname,"w");
01141                                         
01142                                         if (option_data->outts == YUP)
01143                                                 {
01144                                                         option_data->outwritets = fopen (option_data->outnamets,"w");
01145                                                         
01146                                                 }
01147                                         
01148                                         if ((option_data->outwrite == NULL) || (option_data->outlogfile == NULL) ||\
01149                                             (option_data->outwritets == NULL && option_data->outts == YUP) )
01150                                                 {
01151                                                         printf ("\nCould not open ascii output files for writing\n");
01152                                                         exit (1);
01153                                                 }
01154         
01155                                 }
01156         
01157         /* Write out user variables to Logfile */
01158         write_ud (option_data);                 /* writes user data to a file */
01159 
01160         #ifdef ZDBG
01161                 show_ud(option_data, 1);
01162         #endif
01163 
01164 }

00950 {
00951   char message[THD_MAX_NAME];      /* error message */
00952   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00953   int ierror;                         /* number of errors in editing data */
00954 
00955   
00956   /*----- make an empty copy of input dataset -----*/
00957   new_dset = EDIT_empty_copy( dset_time ) ;
00958   
00959   
00960   ierror = EDIT_dset_items( new_dset ,
00961                             ADN_prefix , filename ,
00962                             ADN_label1 , filename ,
00963                             ADN_self_name , filename ,
00964                             ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 
00965                                                            GEN_FUNC_TYPE ,
00966                             ADN_none ) ;
00967   
00968   if( ierror > 0 )
00969     {
00970       sprintf (message,
00971                "*** %d errors in attempting to create output dataset!\n", 
00972                ierror);
00973       FIM_error (message);
00974     }
00975   
00976   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00977     {
00978       sprintf (message,
00979                "Output dataset file %s already exists "
00980                " -- cannot continue!\a\n",
00981                new_dset->dblk->diskptr->header_name);
00982       FIM_error (message);
00983     }
00984   
00985   /*----- deallocate memory -----*/   
00986   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00987   
00988 }

01002 {
01003 
01004   if ((option_data->bucket_filename != NULL)
01005       && (option_data->input1D_filename == NULL))
01006     check_one_output_file (dset_time, option_data->bucket_filename);
01007   
01008 }

00182 {
00183   printf (
00184         "The program estimates the time delay between each voxel time series    \n"
00185         "in a 3D+time dataset and a reference time series[1][2].                \n"
00186         "The estimated delays are relative to the reference time series.\n"
00187         "For example, a delay of 4 seconds means that the voxel time series \n"
00188         "is delayed by 4 seconds with respectto the reference time series.\n\n"
00189         "                                                                       \n"
00190         "Usage:                                                                 \n"
00191         "3ddelay                                                                 \n"
00192         "-input fname       fname = filename of input 3d+time dataset           \n"
00193         "-ideal_file rname  rname = input ideal time series file name           \n"
00194         "   The length of the reference time series should be equal to           \n"
00195         "     that of the 3d+time data set. \n"
00196         "     The reference time series vector is stored in an ascii file.        \n"
00197         "     The programs assumes that there is one value per line and that all  \n"
00198         "     values in the file are part of the reference vector.                \n"
00199         "     PS: Unlike with 3dfim, and FIM in AFNI, values over 33333 are treated\n"
00200         "     as part of the time series.                                          \n" 
00201         "-fs fs             Sampling frequency in Hz. of data time series (1/TR). \n"  
00202         "-T  Tstim          Stimulus period in seconds. \n"
00203         "                   If the stimulus is not periodic, you can set Tstim to 0.\n"
00204         "[-prefix bucket]   The prefix for the results Brick.\n"   
00205         "                   The first subbrick is for Delay.\n"
00206         "                   The second subbrick is for Covariance, which is an estimate\n" 
00207         "                   of the power in voxel time series at the frequencies present \n"
00208         "                   in the reference time series.\n"
00209         "                   The third subbrick is for the Cross Correlation Coefficients between\n"
00210         "                   FMRI time series and reference time series.\n"
00211         "                   The fourth subbrick contains estimates of the Variance of voxel time series.\n"
00212         "                   The default prefix is the prefix of the input 3D+time brick \n"
00213         "                   with a '.DEL' extension appended to it.\n"
00214         "[-uS/-uD/-uR]      Units for delay estimates. (Seconds/Degrees/Radians)\n"
00215         "                   You can't use Degrees or Radians as units unless \n"
00216         "                   you specify a value for Tstim > 0.\n"
00217         "[-phzwrp]          Delay (or phase) wrap.\n"
00218         "                   This switch maps delays from: \n" 
00219         "                   (Seconds) 0->T/2 to 0->T/2 and T/2->T to -T/2->0\n"  
00220         "                   (Degrees) 0->180 to 0->180 and 180->360 to -180->0\n"   
00221         "                   (Radians) 0->pi to 0->pi and pi->2pi to -pi->0\n" 
00222         "                   You can't use this option unless you specify a \n"
00223         "                   value for Tstim > 0.\n\n"
00224         "[-bias]            Do not correct for the bias in the estimates [1][2]\n" 
00225         "[-mask mname]      mname = filename of 3d mask dataset                 \n"
00226         "                   only voxels with non-zero values in the mask would be \n"
00227         "                   considered.                                           \n"
00228         "[-nfirst fnum]     fnum = number of first dataset image to use in      \n"
00229         "                     the delay estimate. (default = 0)                 \n"
00230         "[-nlast  lnum]     lnum = number of last dataset image to use in       \n"
00231         "                     the delay estimate. (default = last)              \n"
00232         "[-nodsamp ]        Do not correct a voxel's estimated delay by the time \n"
00233         "                   at which the slice containing that voxel was acquired.\n\n"
00234         "[-co CCT]          Cross Correlation Coefficient threshold value.\n"
00235         "                   This is only used to limit the ascii output (see below).\n"
00236    "[-nodtrnd]         Do not remove the linear trend from the data time series.\n"
00237    "                   Only the mean is removed. Regardless of this option, \n"
00238    "                   No detrending is done to the reference time series.\n"
00239         "[-asc [out]]       Write the results to an ascii file for voxels with \n"
00240         "[-ascts [out]]     cross correlation coefficients larger than CCT.\n"
00241         "                   If 'out' is not specified, a default name similar \n"
00242         "                   to the default output prefix is used.\n"
00243         "                   -asc, only files 'out' and 'out.log' are written to disk (see ahead)\n"
00244         "                   -ascts, an additional file, 'out.ts', is written to disk (see ahead)\n"
00245         "                   There a 9 columns in 'out' which hold the following values:\n"
00246         "                    1- Voxel Index (VI) : Each voxel in an AFNI brick has a unique index.\n"
00247         "                          Indices map directly to XYZ coordinates.\n"
00248         "                          See AFNI plugin documentations for more info.\n"
00249         "                    2..4- Voxel coordinates (X Y Z): Those are the voxel slice coordinates.\n"
00250         "                          You can see these coordinates in the upper left side of the \n"
00251         "                          AFNI window.To do so, you must first switch the voxel \n"
00252         "                          coordinate units from mm to slice coordinates. \n"
00253         "                          Define Datamode -> Misc -> Voxel Coords ?\n"
00254         "                          PS: The coords that show up in the graph window\n"
00255         "                              could be different from those in the upper left side \n"
00256         "                              of AFNI's main window.\n"
00257         "                    5- Duff : A value of no interest to you. It is preserved for backward \n"
00258         "                          compatibility.\n"
00259         "                    6- Delay (Del) : The estimated voxel delay.\n"
00260         "                    7- Covariance (Cov) : Covariance estimate.\n"
00261         "                    8- Cross Correlation Coefficient (xCorCoef) : Cross Correlation Coefficient.\n"
00262         "                    9- Variance (VTS) : Variance of voxel's time series.\n\n"
00263         "                   The file 'out' can be used as an input to two plugins:\n"
00264         "                     '4Ddump' and '3D+t Extract'\n\n"
00265         "                   The log file 'out.log' contains all parameter settings used for generating \n"
00266         "                   the output brick. It also holds any warnings generated by the plugin.\n"
00267         "                   Some warnings, such as 'null time series ...' , or \n"
00268         "                   'Could not find zero crossing ...' are harmless. '\n"
00269         "                   I might remove them in future versions.\n\n"
00270         "                   A line (L) in the file 'out.ts' contains the time series of the voxel whose\n"
00271         "                   results are written on line (L) in the file 'out'.\n"
00272         "                   The time series written to 'out.ts' do not contain the ignored samples,\n"
00273         "                   they are detrended and have zero mean.\n\n"
00274         "                                                                      \n"
00275         "Random Comments/Advice:\n"
00276         "   The longer you time series, the better. It is generally recomended that\n"
00277         "   the largest delay be less than N/10, N being the length of the time series.\n"
00278         "   The algorithm does go all the way to N/2.\n\n"
00279         "   If you have/find questions/comments/bugs about the plugin, \n"
00280         "   send me an E-mail: ziad@nih.gov\n\n"
00281         "                          Ziad Saad Dec 8 00.\n\n"
00282         "   [1] : Bendat, J. S. (1985). The Hilbert transform and applications to correlation measurements,\n"
00283         "          Bruel and Kjaer Instruments Inc.\n"
00284         "   [2] : Bendat, J. S. and G. A. Piersol (1986). Random Data analysis and measurement procedures, \n"
00285         "          John Wiley & Sons.\n"
00286    "   Author's publications on delay estimation using the Hilbert Transform:\n"
00287    "   [3] : Saad, Z.S., et al., Analysis and use of FMRI response delays. \n"
00288    "         Hum Brain Mapp, 2001. 13(2): p. 74-93.\n"
00289    "   [4] : Saad, Z.S., E.A. DeYoe, and K.M. Ropella, Estimation of FMRI Response Delays. \n"
00290    "         Neuroimage, 2003. 18(2): p. 494-504.\n\n"   
00291     );
00292 
00293   exit(0);
00294 }

01752 {
01753   float fac=0.0 , top ;
01754   
01755   if( MRI_IS_INT_TYPE(otype) ){
01756     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01757     if (top == 0.0)  fac = 0.0;
01758     else  fac = MRI_TYPE_maxval[otype]/top;
01759   }
01760   
01761   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01762   return ( fac );
01763 }

02106         {
02107                 int xpos,ypos,zpos;
02108                 indexTOxyz (option_data, ncall, &xpos , &ypos , &zpos); 
02109 
02110                 switch (option_data->errcode)
02111                         {
02112                                 case ERROR_NOTHINGTODO:
02113                                         fprintf (option_data->outlogfile,"Nothing to do hilbertdelay_V2 call ");
02114                                         break;
02115                                 case ERROR_LARGENSEG:
02116                                         fprintf (option_data->outlogfile,"Number of segments Too Large ");
02117                                         break;
02118                                 case ERROR_LONGDELAY:
02119                                         fprintf (option_data->outlogfile,"Could not find zero crossing before maxdel limit ");
02120                                         break;
02121                                 case ERROR_SERIESLENGTH:
02122                                         fprintf (option_data->outlogfile,"Vectors have different length ");
02123                                         break;
02124                                 case ERROR_NULLTIMESERIES:
02125                                         fprintf (option_data->outlogfile,"Null time series vector ");
02126                                         break;
02127                                 default:
02128                                         fprintf (option_data->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",option_data->errcode);
02129                                         break;
02130                         }       
02131                 fprintf (option_data->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos  );
02132                 return;
02133         }

01267 {
01268   MRI_IMAGE * im;          /* intermediate float data */
01269   float * ar;              /* pointer to float data */
01270   int ts_length;           /* length of input 3d+time data set */
01271   int it;                  /* time index */
01272 
01273 
01274   /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/
01275   im = THD_extract_series (iv, dset_time, 0);
01276 
01277 
01278   /*----- Verify extraction -----*/
01279   if (im == NULL)  FIM_error ("Unable to extract data from 3d+time dataset");
01280 
01281 
01282   /*----- Now extract time series from MRI_IMAGE -----*/
01283   ts_length = DSET_NUM_TIMES (dset_time);
01284   ar = MRI_FLOAT_PTR (im);
01285   for (it = 0;  it < ts_length;  it++)
01286     {
01287       ts_array[it] = ar[it];
01288     }
01289 
01290 
01291   /*----- Release memory -----*/
01292   mri_free (im);   im = NULL;
01293 
01294 }

00170 {
00171   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00172   exit(1);
00173 }

00434 {
00435   int nopt = 1;                     /* input option argument counter */
00436   int ival, index;                  /* integer input */
00437   float fval;                       /* float input */
00438   char message[THD_MAX_NAME];       /* error message */
00439   int k;                            /* ideal time series index */
00440 
00441 
00442   /*----- does user request help menu? -----*/
00443   if (argc < 2 || strcmp(argv[1], "-help") == 0)  display_help_menu();  
00444 
00445   
00446   /*----- initialize the input options -----*/
00447   initialize_options (option_data); 
00448 
00449   
00450   /*----- main loop over input options -----*/
00451   while (nopt < argc )
00452     {
00453 
00454       /*-----   -input filename   -----*/
00455       if (strcmp(argv[nopt], "-input") == 0)
00456                 {
00457                   nopt++;
00458                   if (nopt >= argc)  FIM_error ("need argument after -input ");
00459                   option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME);
00460                   MTEST (option_data->input_filename);
00461                   strcpy (option_data->input_filename, argv[nopt]);
00462                   nopt++;
00463                   continue;
00464                 }
00465 
00466       /*-----   -mask filename   -----*/
00467       if (strcmp(argv[nopt], "-mask") == 0)
00468                 {
00469                   nopt++;
00470                   if (nopt >= argc)  FIM_error ("need argument after -mask ");
00471                   option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME);
00472                   MTEST (option_data->mask_filename);
00473                   strcpy (option_data->mask_filename, argv[nopt]);
00474                   nopt++;
00475                   continue;
00476                 }
00477             
00478 
00479       /*-----   -nfirst num  -----*/
00480       if (strcmp(argv[nopt], "-nfirst") == 0)
00481                 {
00482                   nopt++;
00483                   if (nopt >= argc)  FIM_error ("need argument after -nfirst ");
00484                   sscanf (argv[nopt], "%d", &ival);
00485                   if (ival < 0)
00486                  FIM_error ("illegal argument after -nfirst ");
00487                   option_data->NFirst = ival;
00488                   nopt++;
00489                   continue;
00490                 }
00491 
00492 
00493       /*-----   -nlast num  -----*/
00494       if (strcmp(argv[nopt], "-nlast") == 0)
00495                 {
00496                   nopt++;
00497                   if (nopt >= argc)  FIM_error ("need argument after -nlast ");
00498                   sscanf (argv[nopt], "%d", &ival);
00499                   if (ival < 0)
00500                  FIM_error ("illegal argument after -nlast ");
00501                   option_data->NLast = ival;
00502                   nopt++;
00503                   continue;
00504                 }
00505 
00506       /*-----   -fs num  -----*/
00507       if (strcmp(argv[nopt], "-fs") == 0)
00508                 {
00509                   nopt++;
00510                   if (nopt >= argc)  FIM_error ("need argument after -fs ");
00511                   sscanf (argv[nopt], "%f", &fval);
00512                   if (fval < 0)
00513                  FIM_error ("illegal argument after -fs ");
00514                   option_data->fs = fval;
00515                   nopt++;
00516                   continue;
00517                 }
00518 
00519       /*-----   -T num  -----*/
00520       if (strcmp(argv[nopt], "-T") == 0)
00521                 {
00522                   nopt++;
00523                   if (nopt >= argc)  FIM_error ("need argument after -T ");
00524                   sscanf (argv[nopt], "%f", &fval);
00525                   if (fval < 0)
00526                  FIM_error ("illegal argument after -T ");
00527                   option_data->T = fval;
00528                   nopt++;
00529                   continue;
00530                 }
00531 
00532       /*-----   -ideal_file rname   -----*/
00533       if (strcmp(argv[nopt], "-ideal_file") == 0)
00534                 {
00535                   nopt++;
00536                   if (nopt >= argc)  FIM_error ("need argument after -ideal_file");
00537 
00538                   k = option_data->num_ideal_files;
00539                   if (k+1 > MAX_FILES)
00540                  {
00541                 sprintf (message, "Too many ( > %d ) ideal time series files. ", 
00542                          MAX_FILES);
00543                 FIM_error (message);
00544                  }
00545 
00546                   option_data->ideal_filename[k] 
00547                  = malloc (sizeof(char)*THD_MAX_NAME);
00548                   MTEST (option_data->ideal_filename[k]);
00549                   strcpy (option_data->ideal_filename[k], argv[nopt]);
00550                   option_data->num_ideal_files++;
00551                   nopt++;
00552                   continue;
00553                 }
00554       
00555 
00556       /*-----   -prefix filename   -----*/
00557       if (strcmp(argv[nopt], "-prefix") == 0)
00558                 {
00559                   nopt++;
00560                   if (nopt >= argc)  FIM_error ("need file prefixname after -bucket ");
00561                   option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
00562                   MTEST (option_data->bucket_filename);
00563                   strcpy (option_data->bucket_filename, argv[nopt]);
00564                   nopt++;
00565                   continue;
00566                 }
00567       
00568                 
00569       /*-----   -uS  -----*/
00570       if (strcmp(argv[nopt], "-uS") == 0)
00571                 {
00572                   option_data->unt = METH_SECONDS;
00573                   nopt++;
00574                   continue;
00575                 }
00576 
00577       /*-----   -uR  -----*/
00578       if (strcmp(argv[nopt], "-uR") == 0)
00579                 {
00580                   option_data->unt = METH_RADIANS;
00581                   nopt++;
00582                   continue;
00583                 }
00584 
00585       /*-----   -uD  -----*/
00586       if (strcmp(argv[nopt], "-uD") == 0)
00587                 {
00588                   option_data->unt = METH_DEGREES;
00589                   nopt++;
00590                   continue;
00591                 }
00592 
00593       /*-----   -phzwrp  -----*/
00594       if (strcmp(argv[nopt], "-phzwrp") == 0)
00595                 {
00596                   option_data->wrp = 1;
00597                   nopt++;
00598                   continue;
00599                 }
00600 
00601       /*-----   -bias  -----*/
00602       if (strcmp(argv[nopt], "-bias") == 0)
00603                 {
00604                   option_data->biasrem = 0;
00605                   nopt++;
00606                   continue;
00607                 }
00608 
00609        /*-----   -nodsamp  -----*/
00610       if (strcmp(argv[nopt], "-nodsamp") == 0)
00611                 {
00612                   option_data->Dsamp = 0;
00613                   nopt++;
00614                   continue;
00615                 }
00616 
00617        /*-----   -nodtrnd  -----*/
00618       if (strcmp(argv[nopt], "-nodtrnd") == 0)
00619                 {
00620                   option_data->dtrnd = 0;
00621                   nopt++;
00622                   continue;
00623                 }
00624      
00625        /*-----   -co num -----*/
00626       if (strcmp(argv[nopt], "-co") == 0)
00627                 {
00628                   nopt++;
00629                   if (nopt >= argc)  FIM_error ("need argument after -co");
00630                   sscanf (argv[nopt], "%f", &fval);
00631                   if (fval < 0)
00632                  FIM_error ("illegal argument after -co ");
00633                   option_data->co = fval;
00634                   nopt++;
00635                   continue;
00636                 }
00637 
00638       /*-----   -asc out   -----*/
00639       if (strcmp(argv[nopt], "-asc") == 0)
00640                 {
00641                   nopt++;
00642                   option_data->out = 1;
00643                   if (nopt >= argc)  { 
00644                         option_data->outname = NULL; 
00645                          option_data->outnamelog = NULL;
00646                         
00647                         continue; }
00648                   if (strncmp(argv[nopt], "-", 1) == 0) {
00649                         option_data->outname = NULL; 
00650                         option_data->outnamelog = NULL;
00651                         continue; }
00652                         
00653                   option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
00654                   option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
00655                   
00656                   MTEST (option_data->outname);
00657                   MTEST (option_data->outnamelog);
00658                   strcpy (option_data->outname, argv[nopt]);
00659                   sprintf (option_data->outnamelog, "%s.log", option_data->outname);
00660         
00661                   nopt++;
00662                   continue;
00663                 }
00664 
00665       /*-----   -ascts out   -----*/
00666       if (strcmp(argv[nopt], "-ascts") == 0)
00667                 {
00668                   nopt++;
00669                   option_data->out = 1;
00670                   option_data->outts = 1;
00671                   if (nopt >= argc)  { 
00672                         option_data->outname = NULL;
00673                          option_data->outnamelog = NULL;
00674                         option_data->outnamets = NULL;
00675                         continue; }
00676                   if (strncmp(argv[nopt], "-", 1) == 0) {
00677                          option_data->outnamelog = NULL;
00678                         option_data->outnamets = NULL;
00679                         option_data->outname = NULL; 
00680                         continue; }
00681                         
00682                   option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
00683                   option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
00684                   option_data->outnamets = malloc (sizeof(char)*(THD_MAX_NAME+3));
00685 
00686                   MTEST (option_data->outname);
00687                  MTEST (option_data->outnamets);
00688                   MTEST (option_data->outnamelog);
00689           
00690                   strcpy (option_data->outname, argv[nopt]);
00691                   sprintf (option_data->outnamets, "%s.ts", option_data->outname);
00692                   sprintf (option_data->outnamelog, "%s.log", option_data->outname);
00693                   
00694                   nopt++;
00695                   continue;
00696                 }
00697 
00698       
00699                 
00700                 
00701                 /*----- unknown command -----*/
00702       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00703       FIM_error (message);
00704       
00705     }
00706 
00707 }

02085         {
02086                 *zpos = (int)ncall / (int)(option_data->nxx*option_data->nyy);
02087                 *ypos = (int)(ncall - *zpos * option_data->nxx * option_data->nyy) / (int)option_data->nxx;
02088                 *xpos = ncall - ( *ypos * option_data->nxx ) - ( *zpos * option_data->nxx * option_data->nyy ) ;
02089                 return;
02090         }

00363 {
00364   int is;                     /* index */
00365 
00366 
00367   /*----- Initialize default values -----*/
00368   option_data->fs = 0;
00369   option_data->co = 0.5;
00370   option_data->T = 0;
00371   option_data->unt = METH_SECONDS;
00372   option_data->wrp = 0;
00373   option_data->Navg = 1;
00374   option_data->Nort = 2;
00375   option_data->Nfit = 2;
00376   option_data->Nseg = 1;
00377   option_data->Pover = 0;
00378   option_data->dtrnd = 1;
00379   option_data->biasrem = 1;
00380   option_data->Dsamp = 1;
00381   option_data->outwrite = NULL;
00382   option_data->outwritets = NULL;
00383   option_data->outlogfile = NULL;
00384   option_data->nxx = 64;
00385   option_data->nyy = 64;
00386   option_data->nzz = 20;
00387   option_data->NFirst = 0;
00388   option_data->NLast  = 32767;
00389   option_data->out = 0;
00390   option_data->outts = 0;
00391 
00392         /* Left over from 3dfim+.c remove inthe future, with care !*/
00393   option_data->polort = 1;
00394   option_data->num_ortts = 0;
00395   option_data->num_idealts = 0;
00396   option_data->N      = 0;
00397   option_data->fim_thr = 0.0999;
00398   option_data->cdisp = -1.0;
00399   option_data->q = 0;
00400   option_data->p = 0;
00401   option_data->num_ort_files = 0;
00402   option_data->num_ideal_files = 0;
00403 
00404 
00405 
00406   /*----- Initialize file names -----*/
00407   option_data->input_filename = NULL;
00408   option_data->mask_filename = NULL;  
00409   option_data->input1D_filename = NULL;
00410   option_data->bucket_filename = NULL;
00411   option_data->outname = NULL;  
00412 
00413   for (is = 0;  is < MAX_FILES;  is++)
00414     {  
00415       option_data->ort_filename[is] = NULL;
00416       option_data->ideal_filename[is] = NULL;
00417     }
00418 
00419 }

01229 {
01230   /*----- Allocate memory -----*/
01231   *option_data = (DELAY_options *) malloc (sizeof(DELAY_options));
01232 
01233    
01234   /*----- Get command line inputs -----*/
01235   get_options (argc, argv, *option_data);
01236 
01237 
01238   /*----- Read input data -----*/
01239   read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length,
01240                    ort_array, ort_list, ideal_array, ideal_list);
01241 
01242 
01243   /*----- Check for valid inputs -----*/
01244   check_for_valid_inputs (*option_data, *dset_time, *mask_dset, 
01245                           *fmri_length, ort_array, ideal_array);
01246   
01247 
01248   /*----- Allocate memory for output volumes -----*/
01249   if ((*option_data)->input1D_filename == NULL)
01250     allocate_memory (*option_data, *dset_time, fim_params_vol);
01251 
01252 }

02209 {
02210   DELAY_options * option_data;              /* fim algorithm options */
02211   THD_3dim_dataset * dset_time = NULL;    /* input 3d+time data set */
02212   THD_3dim_dataset * mask_dset = NULL;    /* input mask data set */
02213   float * fmri_data = NULL;               /* input fMRI time series data */
02214   int fmri_length;                        /* length of fMRI time series */
02215   MRI_IMAGE * ort_array[MAX_FILES];       /* ideal time series arrays */
02216   int * ort_list[MAX_FILES];              /* list of ideal time series */
02217   MRI_IMAGE * ideal_array[MAX_FILES];     /* ideal time series arrays */
02218   int * ideal_list[MAX_FILES];            /* list of ideal time series */
02219 
02220   float ** fim_params_vol = NULL;
02221                                 /* array of volumes of fim output parameters */
02222 
02223   
02224   /*----- Identify software -----*/
02225   printf ("\n\n");
02226   printf ("Program: %s \n", PROGRAM_NAME);
02227   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
02228   printf ("Date:    %s \n", PROGRAM_DATE);
02229   printf ("\n");
02230 
02231   /*----- Program initialization -----*/
02232   initialize_program (argc, argv, &option_data, &dset_time, &mask_dset, 
02233                       &fmri_data, &fmri_length, 
02234                       ort_array, ort_list, ideal_array, ideal_list, 
02235                       &fim_params_vol);
02236 
02237 
02238   /*----- Perform fim analysis -----*/
02239   calculate_results (option_data, dset_time, mask_dset, 
02240                      fmri_data, fmri_length,
02241                      ort_array, ort_list, ideal_array, ideal_list, 
02242                      fim_params_vol);
02243   
02244 
02245   /*----- Deallocate memory for input datasets -----*/   
02246   if (dset_time != NULL)  
02247     { THD_delete_3dim_dataset (dset_time, False);  dset_time = NULL; }
02248   if (mask_dset != NULL)  
02249     { THD_delete_3dim_dataset (mask_dset, False);  mask_dset = NULL; }
02250 
02251 
02252   /*----- Write requested output files -----*/
02253   if (option_data->input1D_filename == NULL)
02254     output_results (argc, argv, option_data, fim_params_vol);
02255 
02256 
02257   /*----- Terminate program -----*/
02258   terminate_program (&option_data, 
02259                      ort_array, ort_list, ideal_array, ideal_list, 
02260                      &fim_params_vol); 
02261 
02262 }

01973 {
01974   int q;                    /* number of parameters in baseline model */
01975   int num_idealts;           /* number of ideal time series */
01976   int ib;                   /* sub-brick index */
01977   int is;                   /* ideal index */
01978   int ts_length;            /* length of impulse reponse function */
01979   int N;                    /* number of usable data points */
01980 
01981 
01982   /*----- Initialize local variables -----*/
01983   q = option_data->polort + 1;
01984   num_idealts = option_data->num_idealts;
01985   N = option_data->N;
01986 
01987 
01988   /*----- Write the bucket dataset -----*/
01989   if (option_data->bucket_filename != NULL)
01990   {
01991          write_bucket_data (argc, argv, option_data, fim_params_vol);
01992         }
01993 
01994 }

00830 {
00831   char message[THD_MAX_NAME];    /* error message */
00832 
00833   int num_ort_files;       /* number of ort time series files */
00834   int num_ideal_files;     /* number of ideal time series files */
00835   int is;                  /* time series file index */
00836   int polort;              /* degree of polynomial for baseline model */
00837   int num_ortts;           /* number of ort time series */
00838   int num_idealts;         /* number of ideal time series */
00839   int q;                   /* number of parameters in the baseline model */
00840   int p;                   /* number of parameters in the baseline model 
00841                                                         plus number of ideals */
00842 
00843 
00844   /*----- Initialize local variables -----*/
00845   polort          = option_data->polort;
00846   num_ort_files   = option_data->num_ort_files;
00847   num_ideal_files = option_data->num_ideal_files;
00848 
00849 
00850   /*----- Read the input fMRI measurement data -----*/
00851   if (option_data->input1D_filename != NULL)
00852     {
00853       /*----- Read the input fMRI 1D time series -----*/
00854       *fmri_data = read_one_time_series (option_data->input1D_filename, 
00855                                          fmri_length);
00856       if (*fmri_data == NULL)  
00857         { 
00858           sprintf (message,  "Unable to read time series file: %s", 
00859                    option_data->input1D_filename);
00860           FIM_error (message);
00861         }  
00862       *dset_time = NULL;
00863     }
00864 
00865   else if (option_data->input_filename != NULL)
00866     {
00867       /*----- Read the input 3d+time dataset -----*/
00868       *dset_time = THD_open_one_dataset (option_data->input_filename);
00869       if (!ISVALID_3DIM_DATASET(*dset_time))  
00870         { 
00871           sprintf (message,  "Unable to open data file: %s", 
00872                    option_data->input_filename);
00873           FIM_error (message);
00874         }  
00875       THD_load_datablock ((*dset_time)->dblk);
00876 
00877       if (option_data->mask_filename != NULL)
00878         {
00879           /*----- Read the input mask dataset -----*/
00880           *mask_dset = THD_open_dataset (option_data->mask_filename);
00881           if (!ISVALID_3DIM_DATASET(*mask_dset))  
00882             { 
00883               sprintf (message,  "Unable to open mask file: %s", 
00884                        option_data->mask_filename);
00885               FIM_error (message);
00886             }  
00887           THD_load_datablock ((*mask_dset)->dblk);
00888         }
00889     }
00890   else
00891     FIM_error ("Must specify input measurement data");
00892 
00893   
00894   /*----- Read the input ideal time series files -----*/
00895   for (is = 0;  is < num_ideal_files;  is++)
00896     {
00897       ideal_array[is] = read_time_series (option_data->ideal_filename[is], 
00898                                           &(ideal_list[is]));
00899 
00900       if (ideal_array[is] == NULL)
00901         {
00902           sprintf (message,  "Unable to read ideal time series file: %s", 
00903                    option_data->ideal_filename[is]);
00904           FIM_error (message);
00905         }
00906     }
00907 
00908   
00909   /*----- Count number of ort and number of ideal time series -----*/
00910   num_ortts = 0;
00911   for (is = 0;  is < num_ort_files;  is++)
00912     {
00913       if (ort_list[is] == NULL)
00914         num_ortts += ort_array[is]->ny;
00915       else
00916         num_ortts += ort_list[is][0];
00917     }
00918   q = polort + 1 + num_ortts;
00919 
00920   num_idealts = 0;
00921   for (is = 0;  is < num_ideal_files;  is++)
00922     {
00923       if (ideal_list[is] == NULL)
00924         num_idealts += ideal_array[is]->ny;
00925       else
00926         num_idealts += ideal_list[is][0];
00927     }
00928   p = q + num_idealts;
00929 
00930   option_data->num_ortts = num_ortts;
00931   option_data->num_idealts = num_idealts;
00932   option_data->q = q;
00933   option_data->p = p;
00934 
00935  
00936 }

00722 {
00723   char message[THD_MAX_NAME];    /* error message */
00724   char * cpt;                    /* pointer to column suffix */
00725   char filename[THD_MAX_NAME];   /* time series file name w/o column index */
00726   char subv[THD_MAX_NAME];       /* string containing column index */
00727   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00728                               -- used to read 1D ASCII */
00729   float * far;             /* pointer to MRI_IMAGE floating point data */
00730   int nx;                  /* number of time points in time series */
00731   int ny;                  /* number of columns in time series file */
00732   int iy;                  /* time series file column index */
00733   int ipt;                 /* time point index */
00734   float * ts_data;         /* input time series data */
00735 
00736 
00737   /*----- Read the time series file -----*/
00738   flim = mri_read_1D( ts_filename ) ;
00739   if (flim == NULL)
00740     {
00741       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00742       FIM_error (message);
00743     }
00744 
00745   
00746   /*----- Set pointer to data, and set dimensions -----*/
00747   nx = flim->nx;
00748   ny = flim->ny; iy = 0 ;
00749   if( ny > 1 ){
00750     fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
00751   }
00752 
00753   
00754 
00755   /*----- Save the time series data -----*/
00756   *ts_length = nx;
00757   ts_data = (float *) malloc (sizeof(float) * nx);
00758   MTEST (ts_data);
00759   for (ipt = 0;  ipt < nx;  ipt++)
00760     ts_data[ipt] = far[ipt + iy*nx];   
00761   
00762   
00763   mri_free (flim);  flim = NULL;
00764 
00765   return (ts_data);
00766 }

00781 {
00782   char message[THD_MAX_NAME];    /* error message */
00783   char * cpt;                    /* pointer to column suffix */
00784   char filename[THD_MAX_NAME];   /* time series file name w/o column index */
00785   char subv[THD_MAX_NAME];       /* string containing column index */
00786   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00787                               -- used to read 1D ASCII */
00788   float * far;             /* pointer to MRI_IMAGE floating point data */
00789   int nx;                  /* number of time points in time series */
00790   int ny;                  /* number of columns in time series file */
00791 
00792 
00793   /*----- Read the time series file -----*/
00794 
00795   flim = mri_read_1D(ts_filename) ;
00796   if (flim == NULL)
00797     {
00798       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00799       FIM_error (message);
00800     }
00801 
00802   
00803   far = MRI_FLOAT_PTR(flim);
00804   nx = flim->nx;
00805   ny = flim->ny;
00806   *column_list = NULL;  /* mri_read_1D does column selection */
00807 
00808   return (flim);
00809 }

01309 {
01310   int ip;                   /* parameter index */ 
01311 
01312 
01313   /*----- Saved user requested fim parameters -----*/
01314   for (ip = 0;  ip < NBUCKETS;  ip++)
01315     {
01316       if (fim_params_vol[ip] != NULL)
01317         fim_params_vol[ip][iv]  = fim_params[ip];
01318     }
01319 
01320 }

02001         {
02002                 printf ("\n\nUser Data Values at location :%d\n",loc);
02003                 printf ("option_data->rvec= (1st five elements only)");
02004                 disp_vect (option_data->rvec,5);
02005                 printf ("Input Data Set: %s\n", option_data->input_filename);
02006                 printf ("ASCII output file: %s\n", option_data->outname);
02007                 printf ("ASCII output file (log): %s\n", option_data->outnamelog);
02008                 printf ("ASCII output file (ts): %s\n", option_data->outnamets);
02009                 printf ("Mask Data Set: %s\n", option_data->mask_filename);
02010                 printf ("Reference File Name: %s\n", option_data->ideal_filename[0]);
02011                 printf ("Number of voxels in X direction: %d\n", option_data->nxx);
02012                 printf ("Number of voxels in Y direction: %d\n", option_data->nyy);
02013                 printf ("Number of voxels in Z direction: %d\n", option_data->nzz);
02014                 printf ("option_data->fs= %f\n",option_data->fs);
02015                 printf ("option_data->T= %f\n",option_data->T);
02016                 printf ("option_data->unt= %d\n",option_data->unt);
02017                 printf ("option_data->wrp= %d\n",option_data->wrp);
02018                 printf ("option_data->Navg= %d\n",option_data->Navg);
02019                 printf ("option_data->Nort= %d\n",option_data->Nort);
02020                 printf ("option_data->Nfit= %d\n",option_data->Nfit);
02021                 printf ("option_data->Nseg= %d\n",option_data->Nseg);
02022                 printf ("option_data->Pover= %d\n",option_data->Pover);
02023                 printf ("option_data->dtrnd= %d\n",option_data->dtrnd);
02024                 printf ("option_data->biasrem= %d\n",option_data->biasrem);
02025                 printf ("option_data->Dsamp= %d\n",option_data->Dsamp);
02026                 printf ("option_data->co= %f\n",option_data->co);
02027                 printf ("option_data->ln= %d\n",option_data->ln);
02028                 printf ("option_data->errcode= %d\n",option_data->errcode);
02029                 printf ("option_data->out= %d\n",option_data->out);
02030                 printf ("option_data->outts= %d\n",option_data->outts);
02031                 printf ("Hit enter to continue\a\n\n");
02032                 getchar ();
02033                 return;
02034         }

02163 {
02164   int num_idealts;           /* number of ideal time series */
02165   int ip;                   /* parameter index */
02166   int is;                   /* ideal index */
02167 
02168 
02169   /*----- Initialize local variables -----*/
02170   num_idealts = (*option_data)->num_idealts;
02171 
02172 
02173   /*----- Deallocate memory for option data -----*/   
02174   free (*option_data);  *option_data = NULL;
02175 
02176 
02177   /*----- Deallocate memory for ideal time series -----*/
02178   /*
02179   for (is = 0;  is < num_idealts;  is++)
02180     { free (ideal[is]);  ideal[is] = NULL; } 
02181   */
02182 
02183 
02184   /*----- Deallocate space for volume data -----*/
02185   if (*fim_params_vol != NULL)
02186     {
02187       for (ip = 0;  ip < NBUCKETS;  ip++)
02188         {
02189           if ((*fim_params_vol)[ip] != NULL)
02190             { free ((*fim_params_vol)[ip]);   (*fim_params_vol)[ip] = NULL; }
02191         }
02192 
02193       free (*fim_params_vol);   *fim_params_vol  = NULL; 
02194     }
02195          
02196                  
02197 
02198 }

01826 {
01827   THD_3dim_dataset * old_dset = NULL;      /* prototype dataset */
01828   THD_3dim_dataset * new_dset = NULL;      /* output bucket dataset */
01829   char output_prefix[THD_MAX_NAME];     /* prefix name for bucket dataset */
01830   char output_session[THD_MAX_NAME];    /* directory for bucket dataset */
01831   int nbricks;              /* number of sub-bricks in bucket dataset */
01832   short ** bar = NULL;      /* bar[ib] points to data for sub-brick #ib */
01833 
01834   int brick_type;           /* indicates statistical type of sub-brick */
01835   int brick_coef;           /* regression coefficient index for sub-brick */
01836   char brick_label[THD_MAX_NAME]; /* character string label for sub-brick */
01837 
01838   int ierror;               /* number of errors in editing data */
01839   float * volume;           /* volume of floating point data */
01840 
01841   int N;                    /* number of usable data points */
01842   int q;                    /* number of parameters in the ideal model */
01843   int num_idealts;           /* number of ideal time series */
01844   int ip;                   /* parameter index */
01845   int nxyz;                 /* total number of voxels */
01846   int ibrick;               /* sub-brick index */
01847   int nsam; 
01848   int nfit; 
01849   int nort;                 /* degrees of freedom */
01850   char label[THD_MAX_NAME];   /* general label for sub-bricks */
01851 
01852 
01853   /*----- read prototype dataset -----*/
01854   old_dset = THD_open_one_dataset (option_data->input_filename);
01855 
01856     
01857   /*----- Initialize local variables -----*/
01858   nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;  
01859   num_idealts = option_data->num_idealts;
01860   q = option_data->q;
01861   N = option_data->N;
01862 
01863 
01864   /*----- Calculate number of sub-bricks in the bucket -----*/
01865   nbricks = NBUCKETS;
01866 
01867   strcpy (output_prefix, option_data->bucket_filename);
01868   strcpy (output_session, "./");
01869   
01870  
01871   /*----- allocate memory -----*/
01872   bar  = (short **) malloc (sizeof(short *) * nbricks);
01873   MTEST (bar);
01874   
01875 
01876   /*-- make an empty copy of prototype dataset, for eventual output --*/
01877   new_dset = EDIT_empty_copy (old_dset);
01878 
01879 
01880   /*----- Record history of dataset -----*/
01881   tross_Copy_History( old_dset , new_dset ) ;
01882   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01883   sprintf (label, "Output prefix: %s", output_prefix);
01884   tross_Append_History ( new_dset, label);
01885 
01886   
01887   /*----- delete prototype dataset -----*/ 
01888   THD_delete_3dim_dataset( old_dset , False );  old_dset = NULL ;
01889 
01890 
01891   /*----- Modify some structural properties.  Note that the nbricks
01892           just make empty sub-bricks, without any data attached. -----*/
01893   ierror = EDIT_dset_items (new_dset,
01894                             ADN_prefix,          output_prefix,
01895                             ADN_directory_name,  output_session,
01896                             ADN_type,            HEAD_FUNC_TYPE,
01897                             ADN_func_type,       FUNC_BUCK_TYPE,
01898                             ADN_datum_all,       MRI_short ,   
01899                             ADN_ntt,             0,               /* no time */
01900                             ADN_nvals,           nbricks,
01901                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
01902                             ADN_none ) ;
01903   
01904   if( ierror > 0 )
01905     {
01906       fprintf(stderr, 
01907               "*** %d errors in attempting to create bucket dataset!\n", 
01908               ierror);
01909       exit(1);
01910     }
01911   
01912   if (THD_is_file(DSET_HEADNAME(new_dset))) 
01913     {
01914       fprintf(stderr,
01915               "*** Output dataset file %s already exists--cannot continue!\n",
01916               DSET_HEADNAME(new_dset));
01917       exit(1);
01918     }
01919   
01920 
01921   /*----- Attach individual sub-bricks to the bucket dataset -----*/
01922   ibrick = 0;
01923   for (ip = 0;  ip < NBUCKETS;  ip++)        
01924     {                                 
01925       strcpy (brick_label, DELAY_OUTPUT_TYPE_name[ip]);
01926 
01927       if (ip == COFINDX)
01928                         {
01929                           brick_type = FUNC_COR_TYPE;
01930                           nsam = option_data->ln;  nort = option_data->Nort;
01931                   nfit = option_data->Nfit;
01932                         }
01933                 else
01934                         {
01935                           brick_type = FUNC_FIM_TYPE;
01936                           nsam = 0; nfit = 0; nort = 0;
01937                         }
01938 
01939       volume = fim_params_vol[ip];                
01940       attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01941                         brick_type, brick_label, nsam, nfit, nort, bar);  
01942 
01943       ibrick++;
01944     }
01945 
01946 
01947   /*----- write bucket data set -----*/
01948   THD_load_statistics (new_dset);
01949   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01950   fprintf(stderr,"Wrote bucket dataset ");
01951   fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01952 
01953   
01954   /*----- deallocate memory -----*/   
01955   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01956 
01957 }

02040         {
02041                 fprintf (option_data->outlogfile,"\nLogfile output by Hilbert Delay98 plugin\n");
02042                 fprintf (option_data->outlogfile,"\n\nUser Data Values \n");
02043 
02044                 /* check for NULL filenames          22 July 2005 [rickr] */
02045                 fprintf (option_data->outlogfile,"Input Data Set: %s\n",
02046                          CHECK_NULL_STR(option_data->input_filename));
02047                 fprintf (option_data->outlogfile,"Mask Data Set: %s\n",
02048                          CHECK_NULL_STR(option_data->mask_filename));
02049                 fprintf (option_data->outlogfile,"Ascii output file name: %s\n",
02050                          CHECK_NULL_STR(option_data->outname));
02051                 fprintf (option_data->outlogfile,"Reference File Name: %s\n",
02052                          CHECK_NULL_STR(option_data->ideal_filename[0]));
02053 
02054                 fprintf (option_data->outlogfile,"Number of voxels in X direction: %d\n", option_data->nxx);
02055                 fprintf (option_data->outlogfile,"Number of voxels in Y direction: %d\n", option_data->nyy);
02056                 fprintf (option_data->outlogfile,"Number of voxels in Z direction: %d\n", option_data->nzz);
02057                 fprintf (option_data->outlogfile,"Sampling Frequency = %f\n",option_data->fs);
02058                 fprintf (option_data->outlogfile,"Stimulus Period = %f\n",option_data->T);
02059                 fprintf (option_data->outlogfile,"Delay units = %d\n",option_data->unt);
02060                 fprintf (option_data->outlogfile,"Delay wrap = %d\n",option_data->wrp);
02061                 fprintf (option_data->outlogfile,"Number of segments = %d\n",option_data->Nseg);
02062                 fprintf (option_data->outlogfile,"Length of reference time series = %d\n",option_data->ln);
02063                 fprintf (option_data->outlogfile,"Number of fit parameters = %d\n",option_data->Nfit);
02064                 fprintf (option_data->outlogfile,"Number of nuisance parameters (orts)= %d\n",option_data->Nort);
02065                 fprintf (option_data->outlogfile,"Percent overlap = %d\n",option_data->Pover);
02066                 fprintf (option_data->outlogfile,"Plugin detrending = %d (Always 0, mandatory detrending is performed)\n",option_data->dtrnd);
02067                 fprintf (option_data->outlogfile,"Bias correction = %d\n",option_data->biasrem);
02068                 fprintf (option_data->outlogfile,"Acquisition time correction = %d\n",option_data->Dsamp);
02069                 fprintf (option_data->outlogfile,"Correlation Coefficient Threshold= %f\n",option_data->co);
02070                 fprintf (option_data->outlogfile,"Flag for Ascii output file  = %d\n",option_data->out);
02071                 fprintf (option_data->outlogfile,"Flag for Ascii time series file output = %d\n",option_data->outts);
02072                 fprintf (option_data->outlogfile,"\noption_data->errcode (debugging only)= %d\n\n",option_data->errcode);
02073                 fprintf (option_data->outlogfile,"\nThe format for the output file is the following:\n");
02074                 fprintf (option_data->outlogfile,"VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
02075                 fprintf (option_data->outlogfile,"\nError Log <message> <index> <x> <y> <z>\n\n");
02076                 
02077                 return;
02078         }

02141         {       
02142                 int i;
02143                 
02144                 for (i=0;i<option_data->ln;++i)
02145                         {
02146                                 fprintf (option_data->outwritets, "%f\t",ts[i]);
02147                         }
02148                 fprintf (option_data->outwritets,"\n");
02149         }

02093         {
02094                 *ncall = zpos * ( option_data->nxx*option_data->nyy ) + ypos * option_data->nxx + xpos;
02095                 return;
02096         }

  { "Delay", "Covariance", "Corr. Coef.", "Variance" }