00485 {
00486    hilbert_data_V2 uda,*ud;
00487    MRI_IMAGE * tsim;
00488    MCW_idcode * idc ;                          /* input dataset idcode */
00489    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00490    char *tmpstr , * str , *nprfxstr;                 /* strings from user */
00491    int   ntime, nvec ,nprfx, i;
00492         float * vec , fs , T ;
00493                 
00494         /* Allocate as much character space as Bob specifies in afni.h + a bit more */
00495         
00496         tmpstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
00497         nprfxstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
00498         
00499         if (tmpstr == NULL || nprfxstr == NULL) 
00500                                                                           return "********************\n"
00501                                                                                                 "Could not Allocate\n"
00502                                                                                                 "a teeni weeni bit of\n"
00503                                                                                                 "Memory ! \n"
00504                                                                                                 "********************\n";
00505                                                                                                                 
00506         ud = &uda;              /* ud now points to an allocated space */
00507         ud->errcode = 0;        /*reset error flag */
00508         
00509    /*--------------------------------------------------------------------*/
00510    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00511 
00512    /*--------- go to first input line ---------*/
00513                 
00514    PLUTO_next_option(plint) ;
00515 
00516    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00517    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00518    if( old_dset == NULL )
00519       return "*************************\n"
00520              "Cannot find Input Dataset\n"
00521              "*************************"  ;
00522    
00523    ud->dsetname = DSET_FILECODE (old_dset);
00524         ud->nsamp = DSET_NUM_TIMES (old_dset);
00525         ud->Navg = 1 ;    /* Navg does not play a role for the p value, averaging increases sensitivity */
00526         ud->Nort = PLUTO_get_number(plint) ; /* Should be two by default, for mean and linear trend */
00527         ud->Nfit = 2 ;  /* Always 2 for phase and amplitude for this plugin */
00528         /*--------- go to 2nd input line, input time series ---------*/
00529                 
00530         PLUTO_next_option(plint) ;
00531         
00532         tsim = PLUTO_get_timeseries(plint);
00533         if (tsim == NULL) return "No Timeseries Input";
00534         
00535         ud->ln = (int)tsim -> nx;                                                                       /* number of points in each vector */
00536         nvec    = tsim -> ny;                                                                   /* number of vectors */
00537         ud->rvec   = (float *) MRI_FLOAT_PTR(tsim);     /* vec[i+j*nx] = ith point of jth vector */
00538                                                                                                                                 /* for i=0 .. ntime-1 and j=0 .. nvec-1 */
00539         
00540         if (is_vect_null (ud->rvec,ud->ln) == 1)        /* check if ref vect is all zeroes */
00541                 {
00542                         return "Reference vector is all zeros"; 
00543                 }
00544                 
00545         ud->refname = tsim->name;
00546         ud->ignore = PLUTO_get_number(plint) ;    /* get number item */
00547         
00548         str = PLUTO_get_string(plint) ;
00549    ud->Dsamp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00550    
00551    /*--------- go to 3rd input line, sampling frequency, and stimulus period ---------*/
00552         
00553    PLUTO_next_option(plint) ;
00554    
00555    ud->fs = PLUTO_get_number(plint) ;    /* get number item */
00556    ud->T = PLUTO_get_number(plint) ;    /* get number item */
00557    
00558    ud->co = PLUTO_get_number(plint) ;    /* get number item */
00559    str = PLUTO_get_string(plint) ;
00560    ud->biasrem = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00561    
00562    /*--------- go to 4th input line, delay units and wrp option---------*/
00563                 
00564    PLUTO_next_option(plint) ;
00565 
00566    ud->Nseg = (int)PLUTO_get_number(plint) ;    /* get number item */
00567    ud->Pover = (int)PLUTO_get_number(plint) ;    /* get number item */
00568    
00569    str = PLUTO_get_string(plint) ;                                                      /* get string item (the method) */
00570    ud->unt = (int)PLUTO_string_index( str ,                                     /* find it in list it is from */
00571                                  NUM_METHOD_STRINGS ,
00572                                  method_strings ) ;
00573         
00574         str = PLUTO_get_string(plint) ;  
00575         ud->wrp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00576         
00577    /*--------- go to 5th input line Output prefix ---------*/
00578                 
00579    PLUTO_next_option(plint) ;
00580                 
00581    ud->new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00582         
00583         /* check to see if the field is empty */
00584         if (ud->new_prefix == NULL)
00585                         nprfx = 0;
00586                 else
00587                         nprfx = 1;
00588         /* check if the size is larger than 0. I did not want to check for this unless it's allocated */
00589         if (nprfx == 1 && (int)strlen (ud->new_prefix) == 0)
00590                 nprfx = 0;
00591                 
00592         if (nprfx == 0)         /* now create the new name and make new_prefix point to it */
00593                 {
00594                         sprintf (nprfxstr,"%s.DEL",DSET_PREFIX (old_dset));
00595                         ud->new_prefix = nprfxstr;
00596                         /*printf ("New prefix is set to be : %s\n\a",ud->new_prefix);*/
00597                 }
00598         
00599    if( ! PLUTO_prefix_ok(ud->new_prefix) )      /* check if it is OK */
00600       return "************************\n"
00601              "Output Prefix is illegal\n"
00602              "************************"  ;
00603         
00604         str = PLUTO_get_string(plint) ;                                 /* write delays to file ? */
00605         ud->out = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00606         
00607         ud->strout = PLUTO_get_string(plint) ;                          /* strout is for the outiflename, which will be used after the debugging section */
00608         if (ud->strout == NULL)                                         /* if no output name is given, use the new_prefix */
00609                 {ud->strout = ud->new_prefix;}
00610                 else 
00611                         {       
00612                                 if((int)strlen (ud->strout) == 0) ud->strout = ud->new_prefix;
00613                         }
00614                         
00615         str = PLUTO_get_string(plint) ; 
00616         ud->outts = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00617         
00618         /* ------------------Done with user parameters ---------------------------- */
00619         
00620         ud->nxx = (int)old_dset->daxes->nxx;                            /* get data set dimensions */
00621         ud->nyy = (int)old_dset->daxes->nyy;
00622         ud->nzz = (int)old_dset->daxes->nzz;
00623         
00624         /* No need for users to set these options ...*/
00625         
00626         ud->dtrnd = 0;
00627         
00628         if (ud->ln != (ud->nsamp - ud->ignore))
00629                 {
00630                         ud->errcode = ERROR_BADLENGTH;
00631                         return "***************************\n"
00632                                          "Bad time series length \n"
00633                                          "Check reference time series\n"
00634                                          " or the ignore parameter   \n"
00635                                          "***************************\n";
00636                 }
00637         
00638         if ((ud->unt < 0) || (ud->unt > 2))                                                                             /* unt error Check */
00639                 {
00640          ud->errcode = ERROR_WRONGUNIT;
00641          return "***********************\n"
00642                          " internal error: (ziad)\n"
00643                                          "unt values out of bound\n"
00644                                          "***********************\n";                   /*unt must be between 0 and 2 */
00645                 }
00646           
00647           if ((ud->wrp < 0) || (ud->wrp > 1))                                                                           /* wrp error Check */
00648                 {
00649          ud->errcode = ERROR_WARPVALUES;
00650          return "***********************\n"
00651                          " internal error: (ziad)\n"
00652                                          "wrp values out of bound\n"
00653                                          "***********************\n";                   /* wrp must be between 0 and 1*/
00654                 }
00655           
00656           if (ud->fs < 0.0) {                                         /* fs error Check */
00657          ud->errcode = ERROR_FSVALUES;
00658          return "***********************\n"
00659                          " internal error: (ziad)\n"
00660                                          "fs value is negative !\n"
00661                                          "***********************\n";                   /* fs must be >= 0*/
00662         }
00663           
00664           if (ud->T < 0.0) {                                          /* T error Check */
00665          ud->errcode = ERROR_TVALUES;
00666          return "***********************\n"
00667                          " internal error: (ziad)\n"
00668                                          "T value is negative !\n"
00669                                          "***********************\n";                                   /*T must be >= 0  */
00670         }
00671         
00672                 
00673      if ((ud->T == 0.0) && (ud->unt > 0))                                 /* unt error Check */
00674         {
00675          ud->errcode = ERROR_TaUNITVALUES;
00676          return "***********************\n"
00677                          " internal error: (ziad)\n"
00678                                          "T and unt val. mismatch\n"
00679                                          "***********************\n";                   /*T must be specified, and > 0 in order to use polar units*/
00680         }
00681 
00682     
00683     if ((ud->wrp == 1) && (ud->T == 0.0))                                   /* wrp error Check */
00684         {
00685          ud->errcode = ERROR_TaWRAPVALUES;
00686          return "***********************\n"
00687                          " internal error: (ziad)\n"
00688                                          "wrp and T val. mismatch\n"
00689                                          "***********************\n";                   /*T must be specified, and > 0 in order to use polar warp*/
00690         }
00691          if ((ud->out == NOPE) && (ud->outts == YUP))
00692                         {
00693                          ud->errcode = ERROR_OUTCONFLICT;
00694                          return"***********************\n"
00695                          "error: \n"
00696                                          "Write flag must be on\n"
00697                                          "to use Write ts\n"
00698                                          "***********************\n";   
00699                         
00700                         }
00701         
00702 
00703         /* Open the logfile, regardless of the ascii output files */
00704         sprintf ( tmpstr , "%s.log" , ud->strout);
00705         ud->outlogfile = fopen (tmpstr,"w");
00706 
00707 
00708         if (ud->out == YUP)                                                                     /* open outfile */
00709                                 {                                       
00710                                         ud->outwrite = fopen (ud->strout,"w");
00711                                         
00712                                         if (ud->outts == YUP)
00713                                                 {
00714                                                         sprintf ( tmpstr , "%s.ts" , ud->strout);
00715                                                         ud->outwritets = fopen (tmpstr,"w");
00716                                                         
00717                                                 }
00718                                         
00719                                         if ((ud->outwrite == NULL) || (ud->outlogfile == NULL) ||\
00720                                             (ud->outwritets == NULL && ud->outts == YUP) )
00721                                                 {
00722                                                         ud->errcode = ERROR_FILEOPEN; 
00723                                                         
00724                                                         return "***********************\n"
00725                                                                          "Could Not Write Outfile\n"
00726                                                                          "***********************\n";
00727                                                 }
00728         
00729                                 }
00730         
00731         /* Write out user variables to Logfile */
00732         write_ud (ud);                  /* writes user data to a file */
00733         
00734         /*show_ud (ud,0);       */                      /* For some debugging */
00735 
00736    
00737    /*------------- ready to compute new dataset -----------*/
00738 
00739    new_dset = MAKER_4D_to_typed_fbuc ( old_dset ,             /* input dataset */
00740                                ud->new_prefix ,           /* output prefix */
00741                                -1,                                                      /* negative value indicating data type is like original brick */
00742                                ud->ignore ,               /* ignore count */
00743                                1 ,                    /* detrend = ON Let BOB do it*/
00744                                NBUCKETS,                                        /*Number of values at each voxel*/
00745                                                                                  DELAY_tsfuncV2 ,         /* timeseries processor (bucket version)*/
00746                                                                                  (void *)ud          /* data for tsfunc */
00747                                                                                 ) ; 
00748                                                                                  
00749    /* Setup the label, keywords and types of subbricks */
00750         i = 0;
00751         while (i < NBUCKETS)
00752                 {
00753                         switch (i)
00754                                 {
00755                                         case DELINDX:                                   /* delay value in results vector */
00756                                                 EDIT_BRICK_LABEL (new_dset,i,"Delay");
00757                                                 EDIT_BRICK_ADDKEY (new_dset,i,"D");
00758                                                 ++i;
00759                                                 break;
00760                                         case COVINDX:                                   /* covariance value in results vector */
00761                                                 EDIT_BRICK_LABEL (new_dset,i,"Covariance");
00762                                                 EDIT_BRICK_ADDKEY (new_dset,i,"I");
00763                                                 ++i;
00764                                                 break;
00765                                         case COFINDX:                                   /* cross correlation coefficient value in results vector */
00766                                                 EDIT_BRICK_LABEL (new_dset,i,"Corr. Coef.");
00767                                                 EDIT_BRICK_ADDKEY (new_dset,i,"r");
00768                                                 /* Here you must modify either ud->Nfit or ud->Nort or most likely ud->nsamp based on ud->Navg */
00769                                                 EDIT_BRICK_TO_FICO (new_dset,i,ud->nsamp - ud->ignore,ud->Nfit,ud->Nort);
00770                                                 ++i;
00771                                                 break;
00772                                         case VARINDX:                                   /* FMRI time course variance value in results vector */
00773                                                 EDIT_BRICK_LABEL (new_dset,i,"Variance");
00774                                                 EDIT_BRICK_ADDKEY (new_dset,i,"S2");
00775                                                 ++i;
00776                                                 break;
00777                                         default :
00778                                                 return "*********************\n"
00779                                                                  "Internal Error (ziad)\n"
00780                                                                  " Bad i value \n"
00781                                                                  "*********************\n";
00782                                                 break;
00783                                 }
00784                                 
00785                 }
00786         
00787         PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00788         
00789         
00790 
00791    if (ud->out == YUP)                                                                  /* close outfile and outlogfile*/
00792                                 {
00793                                         fclose (ud->outlogfile);
00794                                         fclose (ud->outwrite);
00795                                         if (ud->outts  == YUP) fclose (ud->outwritets);
00796                                 }
00797                                 else
00798                                 {
00799                                         if (ud->outlogfile != NULL)     fclose (ud->outlogfile);                /* close outlogfile */
00800                                 }
00801         
00802         free (tmpstr);          
00803         free (nprfxstr);
00804    return NULL ;  /* null string returned means all was OK */
00805 }

00814 {
00815    static int nvox , ncall ;
00816         hilbert_data_V2 uda,*ud;
00817         float del, xcorCoef, buckara[4];
00818         float xcor=0.0 ,  tmp=0.0 , tmp2 = 0.0 ,  dtx = 0.0 ,\
00819                          delu = 0.0 , slp = 0.0 , vts = 0.0 , vrvec = 0.0 ;
00820         int i , is_ts_null , status , opt , actv , zpos , ypos , xpos ;
00821         
00822         ud = &uda;
00823         ud = (hilbert_data_V2 *) udp;
00824         
00825    /** is this a "notification"? **/
00826 
00827    if( buckar == NULL ){
00828 
00829       if( npts > 0 ){  /* the "start notification" */
00830 
00831          PLUTO_popup_meter( global_plint ) ;  /* progress meter  */
00832          nvox  = npts ;                       /* keep track of   */
00833          ncall = 0 ;                          /* number of calls */
00834                         
00835       } else {  /* the "end notification" */
00836                         
00837                         opt = 0;                                        /* cleanup in hdelay */
00838                 status = hilbertdelay_V2 (ts,ud->rvec,ud->ln,ud->Nseg,ud->Pover,opt,ud->dtrnd,dtx,ud->biasrem,&tmp,&slp,&xcor,&tmp2,&vts,&vrvec);                                       /* cleanup time */
00839         
00840          PLUTO_set_meter( global_plint , 100 ) ; /* set meter to 100% */
00841 
00842       }
00843       return ;
00844    }
00845 
00846         /* In the old version, I had to check for a match between the lengths of the reference time series and FMRI time series
00847         This is now done before the function is called. */
00848    
00849    if (is_vect_null (ts,npts) == 1) /* check for null vectors */
00850         {
00851                 ud->errcode = ERROR_NULLTIMESERIES;
00852                         error_report (ud , ncall );     /* report the error */
00853                 
00854                 del = 0.0;                                                              /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
00855                 xcorCoef = 0.0;                                         /*  because the data might still be OK */
00856                 xcor = 0.0;
00857         }
00858    
00859    if (ud->errcode == 0)                                /* if there are no errors, proceed */   
00860                 {/* ud->errcode == 0 outer loop */
00861                         opt = 1;                                        /* activate hdelay */
00862                 
00863                 /* transform dtx from seconds to sampling units and correct for the number of points ignored*/
00864                 if (ud->Dsamp == YUP) 
00865                         dtx = (float) (T0 / TR) - ud->ignore;
00866                 else
00867                         dtx = 0.0;
00868                         
00869                 ud->errcode = hilbertdelay_V2 (ts,ud->rvec,ud->ln,ud->Nseg,ud->Pover,opt,ud->dtrnd,dtx,ud->biasrem,&delu,&slp,&xcor,&xcorCoef,&vts,&vrvec);                                     /* cleanup time */
00870                         
00871                         if (ud->errcode == 0) /* If there are no errors, proceed */
00872                                 { /*ud->errcode == 0 inner loop */
00873                                         hunwrap (delu, (float)(1/TR), ud->T, slp, ud->wrp, ud->unt, &del );
00874                                         
00875                                         actv = 1;                                               /* assume voxel is active */
00876         
00877                                         if (xcorCoef < ud->co) actv = 0;                        /* determine if voxel is activated using xcorCoef  */
00878         
00879                                         if ((actv == 1) && (ud->out == YUP))            /* if voxel is truly activated, write results to file without modifying return value */
00880                                                 {
00881                                                         indexTOxyz ( ud , ncall, &xpos , &ypos , &zpos);
00882                                                         fprintf (ud->outwrite,"%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\n", ncall , xpos , ypos , zpos ,  delu , del , xcor , xcorCoef , vts);
00883                                                         if (ud->outts == YUP)
00884                                                                 {
00885                                                                         writets (ud,ts);        
00886                                                                 }
00887                                                 }
00888 
00889                                 }/*ud->errcode == 0 inner loop */
00890                         
00891                         else if (ud->errcode == ERROR_LONGDELAY)
00892                                                 {
00893                                                         error_report ( ud , ncall);     
00894                                                         
00895                                                         del = 0.0;                                                              /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
00896                                                 xcorCoef = 0.0;                                         /*  because the data might still be OK */
00897                                                 xcor = 0.0;
00898                                                         
00899                                                 }
00900                                         else if (ud->errcode != 0)
00901                                                                 {
00902                                                                         error_report ( ud , ncall);     
00903                                                                         
00904                                                                         del = 0.0;                                                              /* Set all the variables to Null and set xcorCoef to an impossible value*/
00905                                                                 xcorCoef = NOWAYXCORCOEF;                                               
00906                                                                 xcor = 0.0;
00907                                                                 }
00908 
00909    }/* ud->errcode == 0 outer loop */
00910    
00911         /* Now fill up the bucket array */
00912         
00913         buckar[DELINDX] = del;
00914         buckar[COVINDX] = xcor;
00915         buckar[COFINDX] = xcorCoef;
00916         buckar[VARINDX] = vts;
00917 
00918 
00919    /** set the progress meter to the % of completion **/
00920    ncall++ ;
00921    
00922    PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
00923    
00924    ud->errcode = 0;     /* Rest error to nothing */
00925    
00926    return ;
00927 }

01029         {
01030                 int xpos,ypos,zpos;
01031                 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos); 
01032 
01033                 switch (ud->errcode)
01034                         {
01035                                 case ERROR_NOTHINGTODO:
01036                                         fprintf (ud->outlogfile,"Nothing to do hilbertdelay_V2 call ");
01037                                         break;
01038                                 case ERROR_LARGENSEG:
01039                                         fprintf (ud->outlogfile,"Number of segments Too Large ");
01040                                         break;
01041                                 case ERROR_LONGDELAY:
01042                                         fprintf (ud->outlogfile,"Could not find zero crossing before maxdel limit ");
01043                                         break;
01044                                 case ERROR_SERIESLENGTH:
01045                                         fprintf (ud->outlogfile,"Vectors have different length ");
01046                                         break;
01047                                 case ERROR_NULLTIMESERIES:
01048                                         fprintf (ud->outlogfile,"Null time series vector ");
01049                                         break;
01050                                 default:
01051                                         fprintf (ud->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",ud->errcode);
01052                                         break;
01053                         }       
01054                 fprintf (ud->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos  );
01055                 return;
01056         }

00353 {
00354     static COMPLEX *w;           /* used to store the w complex array */
00355     static int mstore = 0;       /* stores m for future reference */
00356     static int n = 1;            /* length of fft stored for future */
00357 
00358     COMPLEX u,temp,tm;
00359     COMPLEX *xi,*xip,*xj,*wptr;
00360 
00361     int i,j,k,l,le,windex;
00362 
00363     double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
00364 
00365     if(m != mstore) {
00366 
00367 /* free previously allocated storage and set new m */
00368 
00369         if(mstore != 0) free(w);
00370         mstore = m;
00371         if(m == 0) return;       /* if m=0 then done */
00372 
00373 /* n = 2**m = fft length */
00374 
00375         n = 1 << m;
00376         le = n/2;
00377 
00378 /* allocate the storage for w */
00379 
00380         w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
00381         if(!w) {
00382             printf("\nUnable to allocate complex W array\n");
00383             exit(1);
00384         }
00385 
00386 /* calculate the w values recursively */
00387 
00388         arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
00389         wrecur_real = w_real = cos(arg);
00390         wrecur_imag = w_imag = -sin(arg);
00391         xj = w;
00392         for (j = 1 ; j < le ; j++) {
00393             xj->real = (float)wrecur_real;
00394             xj->imag = (float)wrecur_imag;
00395             xj++;
00396             wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
00397             wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
00398             wrecur_real = wtemp_real;
00399         }
00400     }
00401 
00402 /* start fft */
00403 
00404     le = n;
00405     windex = 1;
00406     for (l = 0 ; l < m ; l++) {
00407         le = le/2;
00408 
00409 /* first iteration with no multiplies */
00410 
00411         for(i = 0 ; i < n ; i = i + 2*le) {
00412             xi = x + i;
00413             xip = xi + le;
00414             temp.real = xi->real + xip->real;
00415             temp.imag = xi->imag + xip->imag;
00416             xip->real = xi->real - xip->real;
00417             xip->imag = xi->imag - xip->imag;
00418             *xi = temp;
00419         }
00420 
00421 /* remaining iterations use stored w */
00422 
00423         wptr = w + windex - 1;
00424         for (j = 1 ; j < le ; j++) {
00425             u = *wptr;
00426             for (i = j ; i < n ; i = i + 2*le) {
00427                 xi = x + i;
00428                 xip = xi + le;
00429                 temp.real = xi->real + xip->real;
00430                 temp.imag = xi->imag + xip->imag;
00431                 tm.real = xi->real - xip->real;
00432                 tm.imag = xi->imag - xip->imag;             
00433                 xip->real = tm.real*u.real - tm.imag*u.imag;
00434                 xip->imag = tm.real*u.imag + tm.imag*u.real;
00435                 *xi = temp;
00436             }
00437             wptr = wptr + windex;
00438         }
00439         windex = 2*windex;
00440     }            
00441 
00442 /* rearrange data by bit reversing */
00443 
00444     j = 0;
00445     for (i = 1 ; i < (n-1) ; i++) {
00446         k = n/2;
00447         while(k <= j) {
00448             j = j - k;
00449             k = k/2;
00450         }
00451         j = j + k;
00452         if (i < j) {
00453             xi = x + i;
00454             xj = x + j;
00455             temp = *xj;
00456             *xj = *xi;
00457             *xi = temp;
00458         }
00459     }
00460 }

00478 {
00479     static COMPLEX *w;           /* used to store the w complex array */
00480     static int mstore = 0;       /* stores m for future reference */
00481     static int n = 1;            /* length of ifft stored for future */
00482 
00483     COMPLEX u,temp,tm;
00484     COMPLEX *xi,*xip,*xj,*wptr;
00485 
00486     int i,j,k,l,le,windex;
00487 
00488     double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
00489     float scale;
00490 
00491     if(m != mstore) {
00492 
00493 /* free previously allocated storage and set new m */
00494 
00495         if(mstore != 0) free(w);
00496         mstore = m;
00497         if(m == 0) return;       /* if m=0 then done */
00498 
00499 /* n = 2**m = inverse fft length */
00500 
00501         n = 1 << m;
00502         le = n/2;
00503 
00504 /* allocate the storage for w */
00505 
00506         w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
00507         if(!w) {
00508             printf("\nUnable to allocate complex W array\n");
00509             exit(1);
00510         }
00511 
00512 /* calculate the w values recursively */
00513 
00514         arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
00515         wrecur_real = w_real = cos(arg);
00516         wrecur_imag = w_imag = sin(arg);  /* opposite sign from fft */
00517         xj = w;
00518         for (j = 1 ; j < le ; j++) {
00519             xj->real = (float)wrecur_real;
00520             xj->imag = (float)wrecur_imag;
00521             xj++;
00522             wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
00523             wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
00524             wrecur_real = wtemp_real;
00525         }
00526     }
00527 
00528 /* start inverse fft */
00529 
00530     le = n;
00531     windex = 1;
00532     for (l = 0 ; l < m ; l++) {
00533         le = le/2;
00534 
00535 /* first iteration with no multiplies */
00536 
00537         for(i = 0 ; i < n ; i = i + 2*le) {
00538             xi = x + i;
00539             xip = xi + le;
00540             temp.real = xi->real + xip->real;
00541             temp.imag = xi->imag + xip->imag;
00542             xip->real = xi->real - xip->real;
00543             xip->imag = xi->imag - xip->imag;
00544             *xi = temp;
00545         }
00546 
00547 /* remaining iterations use stored w */
00548 
00549         wptr = w + windex - 1;
00550         for (j = 1 ; j < le ; j++) {
00551             u = *wptr;
00552             for (i = j ; i < n ; i = i + 2*le) {
00553                 xi = x + i;
00554                 xip = xi + le;
00555                 temp.real = xi->real + xip->real;
00556                 temp.imag = xi->imag + xip->imag;
00557                 tm.real = xi->real - xip->real;
00558                 tm.imag = xi->imag - xip->imag;             
00559                 xip->real = tm.real*u.real - tm.imag*u.imag;
00560                 xip->imag = tm.real*u.imag + tm.imag*u.real;
00561                 *xi = temp;
00562             }
00563             wptr = wptr + windex;
00564         }
00565         windex = 2*windex;
00566     }            
00567 
00568 /* rearrange data by bit reversing */
00569 
00570     j = 0;
00571     for (i = 1 ; i < (n-1) ; i++) {
00572         k = n/2;
00573         while(k <= j) {
00574             j = j - k;
00575             k = k/2;
00576         }
00577         j = j + k;
00578         if (i < j) {
00579             xi = x + i;
00580             xj = x + j;
00581             temp = *xj;
00582             *xj = *xi;
00583             *xi = temp;
00584         }
00585     }
00586 
00587 /* scale all results by 1/n */
00588     scale = (float)(1.0/n);
00589     for(i = 0 ; i < n ; i++) {
00590         x->real = scale*x->real;
00591         x->imag = scale*x->imag;
00592         x++;
00593     }
00594 }

01014         {
01015                 *zpos = (int)ncall / (int)(ud->nxx*ud->nyy);
01016                 *ypos = (int)(ncall - *zpos * ud->nxx * ud->nyy) / (int)ud->nxx;
01017                 *xpos = ncall - ( *ypos * ud->nxx ) - ( *zpos * ud->nxx * ud->nyy ) ;
01018                 return;
01019         }

00298 {
00299    PLUGIN_interface * plint ;     /* will be the output of this routine */
00300 
00301    if( ncall > 0 ) return NULL ;  /* only one interface */
00302 
00303    /*---------------- set titles and call point ----------------*/
00304 
00305    plint = PLUTO_new_interface( "Hilbert Delay98" ,
00306                                 "Time delay between FMRI and reference time series" ,
00307                                 helpstring ,
00308                                 PLUGIN_CALL_VIA_MENU , DELAY_main  ) ;
00309 
00310    global_plint = plint ;  /* make global copy */
00311 
00312    /*--------- 1st line: Input dataset ---------*/
00313 
00314    PLUTO_add_option( plint ,
00315                      "Data" ,  /* label at left of input line */
00316                      "Data" ,  /* tag to return to plugin */
00317                      TRUE       /* is this mandatory? */
00318                    ) ;
00319 
00320    PLUTO_add_dataset(  plint ,
00321                        "3D+time" ,        /* label next to button   */
00322                        ANAT_ALL_MASK ,    /* take only EPI datasets */
00323                        FUNC_ALL_MASK ,    /*  No fim funcs   */
00324                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00325                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00326                     ) ;
00327                                                  
00328         PLUTO_add_number( plint ,
00329                     "Nort" ,  /* label next to chooser */
00330                     1 ,         /* smallest possible value */
00331                     100 ,        /* largest possible value (inactivated for now)*/
00332                     0 ,         /* decimal shift (none in this case) */
00333                     2 ,         /* default value */
00334                     FALSE       /* allow user to edit value? */
00335                   ) ;
00336         
00337    /*---------- 2nd line: Input time series ----------*/
00338    
00339    PLUTO_add_option( plint ,
00340                      "Ref." ,  /* label at left of input line */
00341                      "Ref." ,  /* tag to return to plugin */
00342                      TRUE       /* is this mandatory? */
00343                    ) ;
00344 
00345    PLUTO_add_timeseries(plint,"Ref. Vect."); 
00346    
00347    PLUTO_add_number( plint ,
00348                     "Ignore" ,  /* label next to chooser */
00349                     0 ,         /* smallest possible value */
00350                     50 ,        /* largest possible value (inactivated for now)*/
00351                     0 ,         /* decimal shift (none in this case) */
00352                     0 ,         /* default value */
00353                     FALSE       /* allow user to edit value? */
00354                   ) ;
00355         
00356         PLUTO_add_string( plint ,
00357                      "Dsamp" ,  /*label next to textfield */
00358                      2,yn_strings,  /*   strings to choose among */
00359                      1          /* Default option */
00360                    ) ; 
00361                    
00362    /*---------- 3rd line: sampling frequency ----------*/
00363 
00364    PLUTO_add_option( plint ,
00365                      "Sig." ,  /* label at left of input line */
00366                      "Sig." ,  /* tag to return to plugin */
00367                      TRUE       /* is this mandatory? */
00368                    ) ;
00369 
00370    PLUTO_add_number( plint ,
00371                     "fs in Hz" ,  /* label next to chooser */
00372                     0 ,         /* smallest possible value */
00373                     2000 ,        /* largest possible value */
00374                     1 ,         /* decimal shift (none in this case) */
00375                     5 ,         /* default value */
00376                     TRUE       /* allow user to edit value? */
00377                   ) ;
00378         
00379         PLUTO_add_number( plint ,
00380                     "Tstim sec" ,  /* label next to chooser */
00381                     0.0 ,         /* smallest possible value */
00382                     500 ,        /* largest possible value */
00383                     0 ,         /* decimal shift (none in this case) */
00384                     40 ,         /* default value */
00385                     TRUE       /* allow user to edit value? */
00386                   ) ;
00387 
00388         PLUTO_add_number( plint ,
00389                     "C-Off" ,  /* label next to chooser */
00390                     -10 ,         /* smallest possible value */
00391                     10 ,        /* largest possible value */
00392                     1 ,         /* decimal shift  */
00393                     5 ,         /* default value */
00394                     TRUE       /* allow user to edit value? */
00395                   ) ;
00396    
00397    
00398    PLUTO_add_string( plint ,
00399                      "No-bias" ,  /*label next to textfield */
00400                      2,yn_strings,  /*   strings to choose among */
00401                      1          /* Default option */
00402                    ) ; 
00403                   
00404 
00405         
00406    /*---------- 4th line: Delay Units ----------*/
00407 
00408    PLUTO_add_option( plint ,
00409                      "Alg." ,  /* label at left of input line */
00410                      "Alg." ,  /* tag to return to plugin */
00411                      TRUE        /* is this mandatory? */
00412                    ) ;
00413 
00414    PLUTO_add_number( plint ,
00415                     "N seg." ,  /* label next to chooser */
00416                     1 ,         /* smallest possible value */
00417                     1 ,        /* largest possible value (turned Off for the moment, supporting code is present)*/
00418                     0 ,         /* decimal shift (none in this case) */
00419                     1 ,         /* default value */
00420                     FALSE       /* allow user to edit value? */
00421                   ) ;
00422         
00423         PLUTO_add_number( plint ,
00424                     "% ovrlp" ,  /* label next to chooser */
00425                     0 ,         /* smallest possible value */
00426                     0 ,        /* largest possible value (not implemented)*/
00427                     0 ,         /* decimal shift (none in this case) */
00428                     0 ,         /* default value */
00429                     FALSE       /* allow user to edit value? */
00430                   ) ;
00431 
00432         
00433    PLUTO_add_string( plint ,
00434                      "Units" ,  /* label next to textfield */
00435                      3,method_strings,    /* strings to choose among */
00436                      0          /* Default option */
00437                    ) ;
00438    
00439    PLUTO_add_string( plint ,
00440                      "Phz Wrp" ,  /* label next to textfield */
00441                      2,yn_strings,    /* strings to choose among */
00442                      0          /* Default option */
00443                    ) ;
00444                   
00445 
00446    /*---------- 5th line: Output dataset ----------*/
00447 
00448    PLUTO_add_option( plint ,
00449                      "Output" ,  /* label at left of input line */
00450                      "Output" ,  /* tag to return to plugin */
00451                      TRUE        /* is this mandatory? */
00452                    ) ;
00453 
00454    PLUTO_add_string( plint ,
00455                      "AFNI Prfx" ,  /* label next to textfield */
00456                      0,NULL ,    /* no fixed strings to choose among */
00457                      19          /* 19 spaces for typing in value */
00458                    ) ;
00459         
00460         PLUTO_add_string( plint ,
00461                      "Write" ,  /* label next to textfield */
00462                      2,yn_strings ,    
00463                      1          
00464                    ) ;
00465                    
00466    PLUTO_add_string( plint , "Filename" , 0 , NULL , 19 ) ;
00467    
00468    PLUTO_add_string( plint ,
00469                      "Write ts" ,  /* label next to textfield */
00470                      2,yn_strings ,    
00471                      1          
00472                    ) ;
00473 
00474    /*--------- done with interface setup ---------*/
00475    return plint ;
00476 }

00611 {
00612     static    COMPLEX  *cf;
00613     static    int      mstore = 0;
00614     int       p,num,k,index;
00615     float     Realsum, Realdif, Imagsum, Imagdif;
00616     double    factor, arg;
00617     COMPLEX   *ck, *xk, *xnk, *cx;
00618 
00619 /* First call the fft routine using the x array but with
00620    half the size of the real fft */
00621 
00622     p = m - 1;
00623     cx = (COMPLEX *) x;
00624     fft(cx,p);
00625 
00626 /* Next create the coefficients for recombination, if required */
00627 
00628     num = 1 << p;    /* num is half the real sequence length.  */
00629 
00630     if (m!=mstore){
00631       if (mstore != 0) free(cf);
00632       cf = (COMPLEX *) calloc(num - 1,sizeof(COMPLEX));
00633       if(!cf){
00634         printf("\nUnable to allocate trig recomb coefficients.");
00635         exit(1);
00636       }
00637 
00638       factor = 4.0*atan(1.0)/num;
00639       for (k = 1; k < num; k++){
00640         arg = factor*k;
00641         cf[k-1].real = (float)cos(arg);
00642         cf[k-1].imag = (float)sin(arg);
00643       }
00644     }  
00645 
00646 /* DC component, no multiplies */
00647     y[0].real = cx[0].real + cx[0].imag;
00648     y[0].imag = 0.0;
00649 
00650 /* other frequencies by trig recombination */
00651     ck = cf;
00652     xk = cx + 1;
00653     xnk = cx + num - 1;
00654     for (k = 1; k < num; k++){
00655       Realsum = ( xk->real + xnk->real ) / 2;
00656       Imagsum = ( xk->imag + xnk->imag ) / 2;
00657       Realdif = ( xk->real - xnk->real ) / 2;
00658       Imagdif = ( xk->imag - xnk->imag ) / 2;
00659 
00660       y[k].real = Realsum + ck->real * Imagsum
00661                           - ck->imag * Realdif ;
00662 
00663       y[k].imag = Imagdif - ck->imag * Imagsum
00664                           - ck->real * Realdif ;
00665       ck++;
00666       xk++;
00667       xnk--;
00668     }
00669 }

00934         {
00935                 printf ("\n\nUser Data Values at location :%d\n",loc);
00936                 printf ("ud->dsetname= %s\n",ud->dsetname);
00937                 printf ("ud->refname= %s\n",ud->refname);
00938                 printf ("ud->rvec= (1st three elements only)");
00939                 disp_vect (ud->rvec,3);
00940                 printf ("ud->nxx= %d\n",ud->nxx);
00941                 printf ("ud->nyy= %d\n",ud->nyy);
00942                 printf ("ud->nzz= %d\n",ud->nzz);
00943                 printf ("ud->fs= %f\n",ud->fs);
00944                 printf ("ud->T= %f\n",ud->T);
00945                 printf ("ud->co= %f\n",ud->co);
00946                 printf ("ud->unt= %d\n",ud->unt);
00947                 printf ("ud->wrp= %d\n",ud->wrp);
00948                 printf ("ud->Navg= %d\n",ud->Navg);
00949                 printf ("ud->Nort= %d\n",ud->Nort);
00950                 printf ("ud->Nfit= %d\n",ud->Nfit);
00951                 printf ("ud->Nseg= %d\n",ud->Nseg);
00952                 printf ("ud->Pover= %d\n",ud->Pover);
00953                 printf ("ud->dtrnd= %d\n",ud->dtrnd);
00954                 printf ("ud->biasrem= %d\n",ud->biasrem);
00955                 printf ("ud->Dsamp= %d\n",ud->Dsamp);
00956                 printf ("ud->ln= %d\n",ud->ln);
00957                 printf ("ud->nsamp= %d\n",ud->nsamp);
00958                 printf ("ud->ignore= %d\n",ud->ignore);
00959                 printf ("ud->errcode= %d\n",ud->errcode);
00960                 printf ("ud->new_prefix= %s\n",ud->new_prefix);
00961                 printf ("ud->out= %d\n",ud->out);
00962                 printf ("ud->strout= %s\n",ud->strout);
00963                 printf ("ud->outts= %d\n",ud->outts);
00964                 printf ("Hit enter to continue\a\n\n");
00965                 getchar ();
00966                 return;
00967         }

00974         {
00975                 fprintf (ud->outlogfile,"\nLogfile output by Hilbert Delay98 plugin\n");
00976                 fprintf (ud->outlogfile,"\n\nUser Data Values \n");
00977                 fprintf (ud->outlogfile,"Input data set = %s\n",ud->dsetname);
00978                 fprintf (ud->outlogfile,"Reference file name = %s\n",ud->refname);
00979                 fprintf (ud->outlogfile,"Number of voxels in X dimension = %d\n",ud->nxx);
00980                 fprintf (ud->outlogfile,"Number of voxels in Y dimension = %d\n",ud->nyy);
00981                 fprintf (ud->outlogfile,"Number of voxels in Z dimension = %d\n",ud->nzz);
00982                 fprintf (ud->outlogfile,"Sampling Frequency = %f\n",ud->fs);
00983                 fprintf (ud->outlogfile,"Stimulus Period = %f\n",ud->T);
00984                 fprintf (ud->outlogfile,"Threshold Cut Off value = %f\n",ud->co);
00985                 fprintf (ud->outlogfile,"Delay units = %d\n",ud->unt);
00986                 fprintf (ud->outlogfile,"Delay wrap = %d\n",ud->wrp);
00987                 fprintf (ud->outlogfile,"Number of segments = %d\n",ud->Nseg);
00988                 fprintf (ud->outlogfile,"Number of samples in time series = %d\n",ud->nsamp);
00989                 fprintf (ud->outlogfile,"Ignore = %d\n",ud->ignore);
00990                 fprintf (ud->outlogfile,"Length of reference time series = %d\n",ud->ln);
00991                 fprintf (ud->outlogfile,"Number of fit parameters = %d\n",ud->Nfit);
00992                 fprintf (ud->outlogfile,"Number of nuisance parameters (orts)= %d\n",ud->Nort);
00993                 fprintf (ud->outlogfile,"Percent overlap = %d\n",ud->Pover);
00994                 fprintf (ud->outlogfile,"Plugin detrending = %d (Always 0, mandatory detrending is performed)\n",ud->dtrnd);
00995                 fprintf (ud->outlogfile,"Bias correction = %d\n",ud->biasrem);
00996                 fprintf (ud->outlogfile,"Acquisition time correction = %d\n",ud->Dsamp);
00997                 fprintf (ud->outlogfile,"Prefix for birck output = %s\n",ud->new_prefix);
00998                 fprintf (ud->outlogfile,"Flag for Ascii output file  = %d\n",ud->out);
00999                 fprintf (ud->outlogfile,"Ascii output file name = %s\n",ud->strout);
01000                 fprintf (ud->outlogfile,"Flag for Ascii time series file output = %d\n",ud->outts);
01001                 fprintf (ud->outlogfile,"\nud->errcode (debugging only)= %d\n\n",ud->errcode);
01002                 fprintf (ud->outlogfile,"\nThe format for the output file is the following:\n");
01003                 fprintf (ud->outlogfile,"VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
01004                 fprintf (ud->outlogfile,"\nError Log <message> <index> <x> <y> <z>\n\n");
01005                 
01006                 return;
01007         }

01064         {       
01065                 int i;
01066                 
01067                 for (i=0;i<ud->ln;++i)
01068                         {
01069                                 fprintf (ud->outwritets, "%f\t",ts[i]);
01070                         }
01071                 fprintf (ud->outwritets,"\n");
01072         }