00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "afni.h"
00008 
00009 #ifndef ALLOW_PLUGINS
00010 #  error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012 
00013 /***********************************************************************
00014   Plugin to compute power spectrum of a 3D+time dataset.
00015   This is a moderately complex example, showing how to deal
00016   with different data types on input and output.
00017 ************************************************************************/
00018 
00019 /*------------- string to 'help' the user -------------*/
00020 
00021 static char helpstring[] =
00022   " Purpose: Compute 'Power Spectrum' of a 3D+time dataset.\n"
00023   " Input items are:\n"
00024   "   Input = 3D+time dataset to analyze\n"
00025   "\n"
00026   "   Output: Prefix = Filename prefix for new dataset\n"
00027   "           Datum  = How to store results\n"
00028   "\n"
00029   "   Ignore Count   = How many points to ignore at start\n"
00030   "   Taper Percent  = Amount of data to taper (Hamming)\n"
00031   "   FFT Length     = Fourier transform size to use [N]\n"
00032   "                    (If N > size of data, data will be zero)\n"
00033   "                    (padded. 'shortest' means to use N just)\n"
00034   "                    (above the length of the time series.  )\n"
00035   "\n"
00036   " The output dataset will be stored in the 3D+time format, with\n"
00037   " the 'time' index actually being frequency.  The frequency grid\n"
00038   " spacing will be 1/(N*dt), where N=FFT length and dt = input\n"
00039   " dataset time spacing.\n"
00040   "\n"
00041   " The method used is the simplest known: squared periodogram.\n"
00042   " A single FFT is done (i.e., each point has DOF=2.)\n"
00043 ;
00044 
00045 /*------------- strings for output format -------------*/
00046 
00047 static char * type_strings[]
00048   = { "as Input" , "Byte" , "Short" , "Float" } ;
00049 
00050 #define NUM_TYPE_STRINGS (sizeof(type_strings)/sizeof(char *))
00051 
00052 /*------------- strings for FFT length -------------*/
00053 
00054 static char * fft_strings[] =
00055 #if 0
00056    { "shortest", "32", "64", "128", "256", "512", "1024", "2048", "4096" } ;
00057 #else
00058    /*                    3*       15*      2**      5*    */
00059    { "shortest", "32" ,  "48" ,   "60" ,   "64" ,   "80" ,
00060                          "96" ,  "120" ,  "128" ,  "160" ,
00061                         "192" ,  "240" ,  "256" ,  "320" ,
00062                         "384" ,  "480" ,  "512" ,  "640" ,
00063                         "768" ,  "960" , "1024" , "1280" ,
00064                        "1536" , "1920" , "2048"            } ;
00065 #endif
00066 
00067 #define NUM_FFT_STRINGS (sizeof(fft_strings)/sizeof(char *))
00068 
00069 /*--------------- prototypes for internal routines ---------------*/
00070 
00071 char * POWER_main( PLUGIN_interface * ) ;  /* the entry point */
00072 
00073 #undef ALLOW_TESTING
00074 #ifdef ALLOW_TESTING
00075 PLUGIN_interface * TEST_init(void) ;
00076 char * TEST_main( PLUGIN_interface * ) ;  /* the entry point */
00077 #endif
00078 
00079 /***********************************************************************
00080    Set up the interface to the user:
00081     1) Create a new interface using "PLUTO_new_interface";
00082 
00083     2) For each line of inputs, create the line with "PLUTO_add_option"
00084          (this line of inputs can be optional or mandatory);
00085 
00086     3) For each item on the line, create the item with
00087         "PLUTO_add_dataset" for a dataset chooser,
00088         "PLUTO_add_string"  for a string chooser,
00089         "PLUTO_add_number"  for a number chooser.
00090 ************************************************************************/
00091 
00092 
00093 DEFINE_PLUGIN_PROTOTYPE
00094 
00095 PLUGIN_interface * PLUGIN_init( int ncall )
00096 {
00097    PLUGIN_interface * plint ;     /* will be the output of this routine */
00098 
00099    if( ncall > 1 ) return NULL ;  /* two interfaces */
00100 
00101 #ifdef ALLOW_TESTING
00102    if( ncall == 1 ) return TEST_init() ;
00103 #else
00104    if( ncall == 1 ) return NULL ;
00105 #endif
00106 
00107    /*---------------- set titles and call point ----------------*/
00108 
00109    plint = PLUTO_new_interface( "Power Spectrum" ,
00110                                 "Power Spectrum of a 3D+time Dataset" ,
00111                                 helpstring ,
00112                                 PLUGIN_CALL_VIA_MENU , POWER_main  ) ;
00113 
00114    PLUTO_add_hint( plint , "Power Spectrum of a 3D+time Dataset" ) ;
00115 
00116    PLUTO_set_sequence( plint , "A:newdset:statistics" ) ;
00117 
00118    /*--------- 1st line: Input dataset ---------*/
00119 
00120    PLUTO_add_option( plint ,
00121                      "Input" ,  /* label at left of input line */
00122                      "Input" ,  /* tag to return to plugin */
00123                      TRUE       /* is this mandatory? */
00124                    ) ;
00125 
00126    PLUTO_add_dataset(  plint ,
00127                        "---->>" ,         /* label next to button   */
00128                        ANAT_ALL_MASK ,    /* take any anat datasets */
00129                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00130                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00131                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00132                     ) ;
00133 
00134    /*---------- 2nd line: Output dataset ----------*/
00135 
00136    PLUTO_add_option( plint ,
00137                      "Output" ,  /* label at left of input line */
00138                      "Output" ,  /* tag to return to plugin */
00139                      TRUE        /* is this mandatory? */
00140                    ) ;
00141 
00142    PLUTO_add_string(   plint ,
00143                        "Prefix" ,  /* label next to textfield */
00144                        0,NULL ,    /* no fixed strings to choose among */
00145                        19          /* 19 spaces for typing in value */
00146                    ) ;
00147 
00148    PLUTO_add_string(   plint ,
00149                        "Datum" ,          /* label next to chooser button */
00150                        NUM_TYPE_STRINGS , /* number of strings to choose among */
00151                        type_strings ,     /* list of strings to choose among */
00152                        0                  /* index of default string */
00153                    ) ;
00154 
00155    /*--------- Other lines: Parameters ---------*/
00156 
00157    PLUTO_add_option( plint , "Ignore" , "Ignore" , TRUE ) ;
00158 
00159    PLUTO_add_number( plint ,
00160                      "Count" ,   /* label next to chooser */
00161                      0 ,         /* smallest possible value */
00162                      999 ,       /* largest possible value */
00163                      0 ,         /* decimal shift (none in this case) */
00164                      4 ,         /* default value */
00165                      TRUE        /* allow user to edit value? */
00166                    ) ;
00167 
00168    PLUTO_add_option( plint , "Taper" , "Taper" , TRUE ) ;
00169 
00170    PLUTO_add_number( plint ,
00171                      "Percent" ,    /* label next to chooser */
00172                      0 ,            /* smallest possible value */
00173                      10 ,           /* largest possible value */
00174                      -1 ,           /* decimal shift (1 right == 0 to 100) */
00175                      0 ,            /* default value (with shift == 0) */
00176                      FALSE          /* allow user to edit value? */
00177                    ) ;
00178 
00179    PLUTO_add_option( plint , "FFT" , "FFT" , TRUE ) ;
00180 
00181    PLUTO_add_string( plint ,
00182                      "Length" ,         /* label next to chooser */
00183                      NUM_FFT_STRINGS ,  /* number of strings to choose among */
00184                      fft_strings ,      /* list of strings to choose among */
00185                      0                  /* index of default string */
00186                    ) ;
00187 
00188    /*--------- done with interface setup ---------*/
00189 
00190    return plint ;
00191 }
00192 
00193 /***************************************************************************
00194   Main routine for this plugin (will be called from AFNI).
00195   If the return string is not NULL, some error transpired, and
00196   AFNI will popup the return string in a message box.
00197 ****************************************************************************/
00198 
00199 /*------------------ macros to return workspace at exit -------------------*/
00200 
00201 #undef  FREEUP
00202 #define FREEUP(x) do{ if((x) != NULL){free((x)); (x)=NULL;} } while(0)
00203 
00204 #define FREE_WORKSPACE                              \
00205   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
00206       FREEUP(fout) ; FREEUP(cxar) ; FREEUP(tar)  ;  \
00207       FREEUP(fxar) ; FREEUP(fyar) ; FREEUP(dtr)  ;  \
00208     } while(0)
00209 
00210 /*-------------------------------------------------------------------------*/
00211 
00212 char * POWER_main( PLUGIN_interface * plint )
00213 {
00214    MCW_idcode * idc ;                          /* input dataset idcode */
00215    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00216    char * new_prefix , * str ;                 /* strings from user */
00217    int   new_datum , ignore , nfft , ninp ,    /* control parameters */
00218          old_datum , nuse , ntaper , ktbot ;
00219    float taper ;
00220 
00221    byte   ** bptr  = NULL ;  /* one of these will be the array of */
00222    short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
00223    float  ** fptr  = NULL ;  /* (depending on input datum type) */
00224 
00225    complex * cxar  = NULL ;  /* will be array of data to FFT */
00226    float   * fxar  = NULL ;  /* array loaded from input dataset */
00227    float   * fyar  = NULL ;  /* array loaded from input dataset */
00228    float  ** fout  = NULL ;  /* will be array of output floats */
00229 
00230    float   * tar   = NULL ;  /* will be array of taper coefficients */
00231    float   * dtr   = NULL ;  /* will be array of detrending coeff */
00232 
00233    float dfreq , pfact , phi , xr,xi , yr,yi ;
00234    float x0,x1 , y0,y1 , d0fac,d1fac ;
00235    int   nfreq , nvox , perc , new_units ;
00236    int   istr , ii,iip , ibot,itop , kk , icx ;       /* temp variables */
00237 
00238    /*--------------------------------------------------------------------*/
00239    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00240 
00241    /*--------- go to first input line ---------*/
00242 
00243    PLUTO_next_option(plint) ;
00244 
00245    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00246    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00247    if( old_dset == NULL )
00248       return "*************************\n"
00249              "Cannot find Input Dataset\n"
00250              "*************************"  ;
00251 
00252    /*--------- go to second input line ---------*/
00253 
00254    PLUTO_next_option(plint) ;
00255 
00256    new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00257    if( ! PLUTO_prefix_ok(new_prefix) )      /* check if it is OK */
00258       return "************************\n"
00259              "Output Prefix is illegal\n"
00260              "************************"  ;
00261 
00262    str  = PLUTO_get_string(plint) ;              /* get string item (the datum type) */
00263    istr = PLUTO_string_index( str ,              /* find it in the list it came from */
00264                               NUM_TYPE_STRINGS ,
00265                               type_strings ) ;
00266    switch( istr ){
00267       default:
00268       case 0:
00269          new_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;  /* use old dataset type */
00270       break ;
00271 
00272       case 1: new_datum = MRI_byte  ; break ;  /* assign type of user's choice */
00273       case 2: new_datum = MRI_short ; break ;
00274       case 3: new_datum = MRI_float ; break ;
00275    }
00276 
00277    /*--------- go to next input lines ---------*/
00278 
00279    PLUTO_next_option(plint) ;                 /* skip to next line */
00280    ignore = PLUTO_get_number(plint) ;         /* get number item (ignore) */
00281 
00282    PLUTO_next_option(plint) ;                 /* skip to next line */
00283    taper  = PLUTO_get_number(plint) * 0.01 ;  /* get number item (taper %) */
00284 
00285    /* compute FFT length to use */
00286 
00287    PLUTO_next_option(plint) ;          /* skip to next line */
00288 
00289    str  = PLUTO_get_string(plint) ;    /* get string item for FFT count */
00290    ninp = DSET_NUM_TIMES(old_dset) ;   /* number of values in input */
00291    nuse = ninp - ignore ;              /* number of values to actually use */
00292 
00293    if( nuse < 4 )
00294       return "*****************************\n"
00295              "Not enough time points to FFT\n"
00296              "*****************************"  ;
00297 
00298    if( strcmp(str,fft_strings[0]) == 0 ){
00299 
00300       /*-- get next larger power-of-2 --*/
00301 #if 0
00302       for( nfft=32 ; nfft < nuse ; nfft *= 2 ) ; /* loop until nfft >= nuse */
00303 #else
00304       nfft = csfft_nextup_one35(nuse) ;
00305 #endif
00306 
00307    } else {
00308       nfft = strtol( str , NULL , 10 ) ;  /* just convert string to integer */
00309    }
00310 
00311    /* if the input FFT length is less than the data length,
00312       tell the user and truncate the amount of data to use */
00313 
00314    if( nfft < nuse ){
00315       char str[256] ;
00316 
00317       sprintf( str , "******************************\n"
00318                      "Warning:\n"
00319                      " Number of points in FFT =%4d\n"
00320                      " is less than available data\n"
00321                      " in time series = %d\n"
00322                      "******************************" ,
00323                nfft , nuse ) ;
00324 
00325       PLUTO_popup_transient( plint , str ) ;
00326 
00327       nuse = nfft ;  /* can't use more data than the FFT length */
00328    }
00329 
00330    /* compute the number of output points and the output grid spacing */
00331 
00332    nfreq = nfft / 2 ;                                 /* # frequencies */
00333    dfreq = 1.0 / (nfft * DSET_TIMESTEP(old_dset) ) ;  /* frequency grid */
00334 
00335    switch( DSET_TIMEUNITS(old_dset) ){
00336       case UNITS_MSEC_TYPE: dfreq *= 1000.0 ; new_units = UNITS_HZ_TYPE ; break;
00337       case UNITS_SEC_TYPE:                    new_units = UNITS_HZ_TYPE ; break;
00338       case UNITS_HZ_TYPE:                     new_units = UNITS_SEC_TYPE; break;
00339 
00340       default: new_units = DSET_TIMEUNITS(old_dset) ; break ; /* shouldn't happen */
00341    }
00342 
00343    /*------------------------------------------------------*/
00344    /*---------- At this point, the inputs are OK ----------*/
00345 
00346    PLUTO_popup_meter( plint ) ;  /* popup a progress meter */
00347 
00348    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00349 
00350    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
00351 
00352    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
00353 
00354    switch( old_datum ){  /* pointer type depends on input datum type */
00355 
00356       default:
00357          return "******************************\n"
00358                 "Illegal datum in Input Dataset\n"
00359                 "******************************"  ;
00360 
00361       /** create array of pointers into old dataset sub-bricks **/
00362       /** Note that we skip the first 'ignore' sub-bricks here **/
00363 
00364       /*--------- input is bytes ----------*/
00365       /* voxel #i at time #k is bptr[k][i] */
00366       /* for i=0..nvox-1 and k=0..nuse-1.  */
00367 
00368       case MRI_byte:
00369          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
00370          if( bptr == NULL ) return "Malloc\nFailure!\n [bptr]" ;
00371          for( kk=0 ; kk < nuse ; kk++ )
00372             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
00373       break ;
00374 
00375       /*--------- input is shorts ---------*/
00376       /* voxel #i at time #k is sptr[k][i] */
00377       /* for i=0..nvox-1 and k=0..nuse-1.  */
00378 
00379       case MRI_short:
00380          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
00381          if( sptr == NULL ) return "Malloc\nFailure!\n [sptr]" ;
00382          for( kk=0 ; kk < nuse ; kk++ )
00383             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
00384       break ;
00385 
00386       /*--------- input is floats ---------*/
00387       /* voxel #i at time #k is fptr[k][i] */
00388       /* for i=0..nvox-1 and k=0..nuse-1.  */
00389 
00390       case MRI_float:
00391          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
00392          if( fptr == NULL ) return "Malloc\nFailure!\n [fptr]" ;
00393          for( kk=0 ; kk < nuse ; kk++ )
00394             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
00395       break ;
00396 
00397    } /* end of switch on input type */
00398 
00399    /*---- allocate space for 2 voxel timeseries and 1 FFT ----*/
00400 
00401    cxar = (complex *) malloc( sizeof(complex) * nfft ) ; /* FFT */
00402    fxar = (float *)   malloc( sizeof(float) * nuse ) ;   /* input */
00403    fyar = (float *)   malloc( sizeof(float) * nuse ) ;   /* input */
00404    if( cxar == NULL || fxar == NULL || fyar == NULL ){
00405       FREE_WORKSPACE ;
00406       return "Malloc\nFailure!\n [cxar]" ;
00407    }
00408 
00409    /*--------- make space for taper coefficient array ---------*/
00410 
00411    tar = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
00412    dtr = (float *) malloc( sizeof(float) * nuse ) ;
00413 
00414    if( tar == NULL || dtr == NULL ){
00415       FREE_WORKSPACE ;
00416       return "Malloc\nFailure!\n [tar]" ;
00417    }
00418 
00419    ntaper = (int)(0.5 * taper * nuse + 0.49) ; /* will taper data over */
00420    phi    = PI / MAX(ntaper,1) ;               /* kk=0..ntaper-1 on left */
00421    ktbot  = nuse - ntaper ;                    /* kk=ktbot..nuse-1 on right */
00422    pfact  = 0.0 ;                              /* sum of taper**2 */
00423 
00424    for( kk=0 ; kk < nuse ; kk++ ){                       /* Hamming-ize */
00425       if( kk < ntaper )
00426          tar[kk] = 0.54 - 0.46 * cos(kk*phi) ;           /* ramp up */
00427       else if( kk >= ktbot )
00428          tar[kk] = 0.54 + 0.46 * cos((kk-ktbot+1)*phi) ; /* ramp down */
00429       else
00430          tar[kk] = 1.0 ;                                 /* in the middle */
00431 
00432       pfact  += tar[kk] * tar[kk] ;
00433 
00434       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* factors for linear detrending */
00435    }
00436 
00437    d0fac = 1.0 / nuse ;
00438    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
00439 
00440    /*--- compute factor to go from |FFT|**2 to PSD;
00441          includes the scaling needed for loss of energy with tapering ---*/
00442 
00443    pfact = DSET_TIMESTEP(old_dset) / pfact ;
00444 
00445    /*--- include scaling factors for sub-bricks, if any ---*/
00446 
00447    for( kk=0 ; kk < nuse ; kk++ )
00448       if( DSET_BRICK_FACTOR(old_dset,kk+ignore) > 0.0 )
00449          tar[kk] *= DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
00450 
00451    /*---------------------- make a new dataset ----------------------*/
00452 
00453    new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00454 
00455    { char * his = PLUTO_commandstring(plint) ;
00456      tross_Copy_History( old_dset , new_dset ) ;
00457      tross_Append_History( new_dset , his ) ; free(his) ;
00458    }
00459 
00460    /*-- edit some of its internal parameters --*/
00461 
00462    ii = EDIT_dset_items(
00463            new_dset ,
00464               ADN_prefix      , new_prefix ,           /* filename prefix */
00465               ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00466               ADN_datum_all   , new_datum ,            /* atomic datum */
00467               ADN_nvals       , nfreq ,                /* # sub-bricks */
00468               ADN_ntt         , nfreq ,                /* # time points */
00469               ADN_ttorg       , dfreq ,                /* time origin */
00470               ADN_ttdel       , dfreq ,                /* time step */
00471               ADN_ttdur       , dfreq ,                /* time duration */
00472               ADN_nsl         , 0 ,                    /* z-axis time slicing */
00473               ADN_tunits      , new_units ,            /* time units */
00474            ADN_none ) ;
00475 
00476    if( ii != 0 ){
00477       THD_delete_3dim_dataset( new_dset , False ) ;
00478       FREE_WORKSPACE ;
00479       return "***********************************\n"
00480              "Error while creating output dataset\n"
00481              "***********************************"  ;
00482    }
00483 
00484    /*------ make floating point output sub-bricks
00485             (only at the end will scale to byte or shorts)
00486 
00487             Output #ii at freq #kk will go into fout[kk][ii],
00488             for kk=0..nfreq-1, and for ii=0..nvox-1.          ------*/
00489 
00490    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00491 
00492    fout = (float **) malloc( sizeof(float *) * nfreq ) ;  /* ptrs to sub-bricks */
00493 
00494    if( fout == NULL ){
00495       THD_delete_3dim_dataset( new_dset , False ) ;
00496       FREE_WORKSPACE ;
00497       return "Malloc\nFailure!\n [fout]" ;
00498    }
00499 
00500    for( kk=0 ; kk < nfreq ; kk++ ){
00501       fout[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
00502       if( fout[kk] == NULL ) break ;
00503    }
00504 
00505    if( kk < nfreq ){
00506       for( ; kk >= 0 ; kk-- ) FREEUP(fout[kk]) ;   /* free all we did get */
00507       THD_delete_3dim_dataset( new_dset , False ) ;
00508       FREE_WORKSPACE ;
00509       return "Malloc\nFailure!\n [arrays]" ;
00510    }
00511 
00512    { char buf[128] ;
00513      ii = (nfreq * nvox * sizeof(float)) / (1024*1024) ;
00514      sprintf( buf , "  \n"
00515                     "*** 3D+time Power Spectral Density:\n"
00516                     "*** Using %d MBytes of workspace,\n "
00517                     "*** with FFT length = %d\n" , ii,nfft ) ;
00518      PLUTO_popup_transient( plint , buf ) ;
00519    }
00520 
00521    /*----------------------------------------------------*/
00522    /*----- Setup has ended.  Now do some real work. -----*/
00523 
00524    /***** loop over voxels *****/
00525 
00526    for( ii=0 ; ii < nvox ; ii += 2 ){  /* 2 time series at a time */
00527 
00528       iip = (ii+1) % nvox ;           /* voxel after ii */
00529 
00530       /*** load data from input dataset, depending on type ***/
00531 
00532       switch( old_datum ){
00533 
00534          /*** input = bytes ***/
00535 
00536          case MRI_byte:
00537             for( kk=0 ; kk < nuse ; kk++ ){
00538                fxar[kk] = bptr[kk][ii] ;
00539                fyar[kk] = bptr[kk][iip] ;
00540             }
00541          break ;
00542 
00543          /*** input = shorts ***/
00544 
00545          case MRI_short:
00546             for( kk=0 ; kk < nuse ; kk++ ){
00547                fxar[kk] = sptr[kk][ii] ;
00548                fyar[kk] = sptr[kk][iip] ;
00549             }
00550          break ;
00551 
00552          /*** input = floats ***/
00553 
00554          case MRI_float:
00555             for( kk=0 ; kk < nuse ; kk++ ){
00556                fxar[kk] = fptr[kk][ii] ;
00557                fyar[kk] = fptr[kk][iip] ;
00558             }
00559          break ;
00560 
00561       } /* end of switch over input type */
00562 
00563       /*** detrend:
00564              x0 = sum( fxar[kk] )
00565              x1 = sum( fxar[kk] * (kk-0.5*(N-1)) )
00566            x0 is used to remove the mean of fxar
00567            x1 is used to remove the linear trend of fxar ***/
00568 
00569       x0 = x1 = y0 = y1 = 0.0 ;
00570       for( kk=0 ; kk < nuse ; kk++ ){
00571          x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
00572          y0 += fyar[kk] ; y1 += fyar[kk] * dtr[kk] ;
00573       }
00574 
00575       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
00576       y0 *= d0fac ; y1 *= d1fac ;
00577 
00578       for( kk=0 ; kk < nuse ; kk++ ){
00579          fxar[kk] -= (x0 + x1 * dtr[kk]) ;  /* remove mean and trend here! */
00580          fyar[kk] -= (y0 + y1 * dtr[kk]) ;
00581       }
00582 
00583       /*** taper, scale, and put into cxar array ***/
00584 
00585       for( kk=0 ; kk < nuse ; kk++ ){
00586          cxar[kk].r = fxar[kk] * tar[kk] ;
00587          cxar[kk].i = fyar[kk] * tar[kk] ;
00588       }
00589 
00590       /*** load zeros after where data was put ***/
00591 
00592       for( kk=nuse ; kk < nfft ; kk++ ) cxar[kk].r = cxar[kk].i = 0.0 ;
00593 
00594       /***** do the FFT (at long last) *****/
00595 
00596       csfft_cox( -1 , nfft , cxar ) ;
00597 
00598       /***** now compute output into corresponding voxels in fout *****/
00599 
00600       /*--- Let x = fxar (1st real time series)
00601                 y = fyar (2nd real time series)
00602                 z = cxar (complex time series) = x + i y
00603                 N = nfft (length of FFT)
00604 
00605             Then after FFT, since x and y are real, we have
00606               zhat[k]  = xhat[k] + i yhat[k]  > for k=1..N/2
00607             zhat[N-k]* = xhat[k] - i yhat[k]
00608 
00609             so we can untangle the FFTs of x and y by
00610               xhat[k] = 0.5 ( zhat[k] + zhat[N-k]* )
00611               yhat[k] = 0.5 ( zhat[k] - zhat[N-k]* ) / i
00612 
00613             This is the basis for doing 2 time series at once. ---*/
00614 
00615       for( kk=1 ; kk <= nfreq ; kk++ ){
00616          xr = 0.5 * ( cxar[kk].r + cxar[nfft-kk].r ) ; /* Re xhat[kk] */
00617          xi = 0.5 * ( cxar[kk].i - cxar[nfft-kk].i ) ; /* Im xhat[kk] */
00618          yr = 0.5 * ( cxar[kk].i + cxar[nfft-kk].i ) ; /* Re yhat[kk] */
00619          yi = 0.5 * ( cxar[kk].r - cxar[nfft-kk].r ) ; /*-Im yhat[kk] */
00620 
00621          fout[kk-1][ii]  = pfact * (xr*xr + xi*xi) ;
00622          fout[kk-1][iip] = pfact * (yr*yr + yi*yi) ;
00623       }
00624 
00625       perc = (100 * ii) / nvox ;        /* display percentage done */
00626       PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00627 
00628    } /* end of outer loop over 2 voxels at a time */
00629 
00630    DSET_unload( old_dset ) ;  /* don't need this no more */
00631 
00632    /*------------------------------------------------------------*/
00633    /*------- The output is now in fout[kk][ii],
00634              for kk=0..nfreq-1 , ii=0..nvox-1.
00635              We must now put this into the output dataset -------*/
00636 
00637    switch( new_datum ){
00638 
00639       /*** output is floats is the simplest:
00640            we just have to attach the fout bricks to the dataset ***/
00641 
00642       case MRI_float:
00643          for( kk=0 ; kk < nfreq ; kk++ )
00644             EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
00645       break ;
00646 
00647       /*** output is shorts:
00648            we have to create a scaled sub-brick from fout ***/
00649 
00650       case MRI_short:{
00651          short * bout ;
00652          float fac ;
00653 
00654          for( kk=0 ; kk < nfreq ; kk++ ){  /* loop over sub-bricks */
00655 
00656             /*-- get output sub-brick --*/
00657 
00658             bout = (short *) malloc( sizeof(short) * nvox ) ;
00659             if( bout == NULL ){
00660                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
00661                return("\nFinal malloc error in plug_power!\n") ;
00662                /* EXIT(1) ;*/
00663             }
00664 
00665             /*-- find scaling and then scale --*/
00666 
00667             fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
00668             if( fac > 0.0 ){
00669                fac = 32767.0 / fac ;
00670                EDIT_coerce_scale_type( nvox,fac ,
00671                                        MRI_float,fout[kk] , MRI_short,bout ) ;
00672                fac = 1.0 / fac ;
00673             }
00674 
00675             free( fout[kk] ) ;  /* don't need this anymore */
00676 
00677             /*-- put output brick into dataset, and store scale factor --*/
00678 
00679             EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
00680             tar[kk] = fac ;
00681 
00682             perc = (100 * kk) / nfreq ;
00683             PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00684          }
00685 
00686          /*-- save scale factor array into dataset --*/
00687 
00688          EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
00689 
00690       }
00691       break ;
00692 
00693       /*** output is bytes (byte = unsigned char)
00694            we have to create a scaled sub-brick from fout ***/
00695 
00696       case MRI_byte:{
00697          byte * bout ;
00698          float fac ;
00699 
00700          for( kk=0 ; kk < nfreq ; kk++ ){  /* loop over sub-bricks */
00701 
00702             /*-- get output sub-brick --*/
00703 
00704             bout = (byte *) malloc( sizeof(byte) * nvox ) ;
00705             if( bout == NULL ){
00706                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
00707                return("\nFinal malloc error in plug_power!\n\a") ;
00708                /* EXIT(1) ;*/
00709             }
00710 
00711             /*-- find scaling and then scale --*/
00712 
00713             fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
00714             if( fac > 0.0 ){
00715                fac = 255.0 / fac ;
00716                EDIT_coerce_scale_type( nvox,fac ,
00717                                        MRI_float,fout[kk] , MRI_byte,bout ) ;
00718                fac = 1.0 / fac ;
00719             }
00720 
00721             free( fout[kk] ) ;  /* don't need this anymore */
00722 
00723             /*-- put output brick into dataset, and store scale factor --*/
00724 
00725             EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
00726             tar[kk] = fac ;
00727 
00728             perc = (100 * kk) / nfreq ;
00729             PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00730          }
00731 
00732          /*-- save scale factor array into dataset --*/
00733 
00734          EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
00735       }
00736       break ;
00737 
00738    } /* end of switch on output data type */
00739 
00740    /*-------------- Cleanup and go home ----------------*/
00741 
00742    PLUTO_set_meter( plint , 100 ) ;  /* set progress meter to 100% */
00743 
00744    PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00745 
00746    FREE_WORKSPACE ;
00747    return NULL ;  /* null string returned means all was OK */
00748 }
00749 
00750 #ifdef ALLOW_TESTING
00751 /*****************************************************************************
00752  -----------------------------------------------------------------------------
00753            Create the second interface within this plugin.
00754  -----------------------------------------------------------------------------*/
00755 
00756 PLUGIN_interface * TEST_init( void )
00757 {
00758    PLUGIN_interface * plint ;     /* will be the output of this routine */
00759 
00760    /*---------------- set titles and call point ----------------*/
00761 
00762    plint = PLUTO_new_interface( "Testing" ,
00763                                 "Testing, Testing, 1-2-3 ..." ,
00764                                 NULL ,
00765                                 PLUGIN_CALL_VIA_MENU , TEST_main  ) ;
00766 
00767    PLUTO_add_hint( plint , "1-2-3, 1-2-3, ..." ) ;
00768 
00769    PLUTO_add_option( plint ,
00770                      "Input" ,  /* label at left of input line */
00771                      "Input" ,  /* tag to return to plugin */
00772                      TRUE       /* is this mandatory? */
00773                    ) ;
00774 
00775    PLUTO_add_dataset_list(  plint ,
00776                             "Datasets" ,       /* label next to button   */
00777                             ANAT_ALL_MASK ,    /* take any anat datasets */
00778                             FUNC_FIM_MASK ,    /* only allow fim funcs   */
00779                             DIMEN_4D_MASK |    /* need 3D+time datasets  */
00780                             BRICK_ALLREAL_MASK /* need real-valued datasets */
00781                          ) ;
00782    return plint ;
00783 }
00784 
00785 char * TEST_main( PLUGIN_interface * plint )
00786 {
00787    MRI_IMAGE * tsim ;
00788    MCW_idclist * idclist ;
00789    MCW_idcode * idc ;
00790    THD_3dim_dataset * dset ;
00791    char str[256] ;
00792    int id ;
00793 
00794    /*--------- go to first input line ---------*/
00795 
00796    PLUTO_next_option(plint) ;
00797 
00798    idclist = PLUTO_get_idclist(plint) ;
00799    if( PLUTO_idclist_count(idclist) == 0 )
00800       return " \nNo input dataset list!\n " ;
00801 
00802    id = 0 ;
00803    do {
00804       idc  = PLUTO_idclist_next(idclist) ;
00805       dset = PLUTO_find_dset(idc) ;
00806       if( dset == NULL ) return NULL ;
00807       id++ ;
00808       sprintf(str, " \nDataset %d = %s\n nx = %d\n ny = %d\n nz = %d\n " ,
00809               id , DSET_FILECODE(dset) , dset->daxes->nxx,dset->daxes->nyy,dset->daxes->nzz ) ;
00810 
00811       PLUTO_popup_transient( plint , str ) ;
00812    } while(1) ;
00813    return NULL ;
00814 }
00815 #endif  /* ALLOW_TESTING */