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 provide a least squares fitting 1D function for graphs
00015 ************************************************************************/
00016 
00017 /*------------- string to 'help' the user -------------*/
00018 
00019 static char helpstring[] =
00020    " Purpose: Control the 'LSqFit' and 'LSqDtr 1D functions.\n"
00021    "\n"
00022    " Parameters:  Baseline = 'Constant' or 'Linear'\n"
00023    "                           Is the baseline 'a' or 'a+b*t'?\n"
00024    "              Ignore   = Number of points to ignore at\n"
00025    "                           start of each timeseries.\n"
00026    " \n"
00027    " Sinusoids:   Period    = Fundamental period to use.\n"
00028    "              Harmonics = Number of overtones to use.\n"
00029    " \n"
00030    " Timeseries:  File      = Input timeseries file to use.\n"
00031 ;
00032 
00033 static char plehstring[] =
00034    " Purpose: Generate a timeseries and store in AFNI list\n"
00035    "\n"
00036    " Parameters:  Delta  = time step between points\n"
00037    "              Length = number of points\n"
00038    "\n"
00039    " Output:      Label  = String to label timeseries by\n"
00040    "                        in AFNI choosers\n"
00041    "\n"
00042    " Polynomial:  Order  = Maximum power to include\n"
00043    "                       (Chebyshev polynomials are used)\n"
00044    "\n"
00045    " Sinusoid:    Period    = Fundamental period to use.\n"
00046    "              Harmonics = Number of overtones to use.\n"
00047 ;
00048 
00049 /*------- Strings for baseline control ------*/
00050 
00051 #define NBASE 3
00052 static char * baseline_strings[NBASE] = { "Constant" , "Linear" , "Quadratic" } ;
00053 
00054 /*--------------- prototypes for internal routines ---------------*/
00055 
00056 char * LSQ_main( PLUGIN_interface * ) ;  /* the entry point */
00057 
00058 void LSQ_fitter( int nt, double to, double dt, float * vec, char ** label ) ;
00059 void LSQ_detrend( int nt, double to, double dt, float * vec, char ** label ) ;
00060 void LSQ_worker( int nt, double dt, float * vec, int dofit, char ** label ) ;
00061 
00062 PLUGIN_interface * TSGEN_init(void) ;
00063 char * TSGEN_main( PLUGIN_interface * ) ;
00064 
00065 PLUGIN_interface * EXP0D_init(void) ;
00066 char * EXP0D_main( PLUGIN_interface * ) ;
00067 void EXP0D_worker( int num , float * vec ) ;
00068 
00069 #ifdef ALLOW_LOMO
00070 PLUGIN_interface * LOMOR_init(void) ;
00071 char * LOMOR_main( PLUGIN_interface * ) ;
00072 void LOMOR_fitter() ;
00073 #endif
00074 
00075 /*---------------- global data -------------------*/
00076 
00077 static PLUGIN_interface * global_plint = NULL ;
00078 
00079 #define NRMAX_SIN 2
00080 #define NRMAX_TS  2
00081 #define HARM_MAX  22
00082 
00083 static int polort=1 , ignore=3 , nrsin=0 , nrts=0 , initialize=1 ;
00084 static float sinper[NRMAX_SIN] ;
00085 static int   sinharm[NRMAX_SIN] ;
00086 static MRI_IMAGE * tsim[NRMAX_TS] ;
00087 
00088 /***********************************************************************
00089    Set up the interface to the user:
00090     1) Create a new interface using "PLUTO_new_interface";
00091 
00092     2) For each line of inputs, create the line with "PLUTO_add_option"
00093          (this line of inputs can be optional or mandatory);
00094 
00095     3) For each item on the line, create the item with
00096         "PLUTO_add_dataset"    for a dataset chooser,
00097         "PLUTO_add_string"     for a string chooser,
00098         "PLUTO_add_number"     for a number chooser,
00099         "PLUTO_add_timeseries" for a timeseries chooser.
00100 ************************************************************************/
00101 
00102 
00103 DEFINE_PLUGIN_PROTOTYPE
00104 
00105 PLUGIN_interface * PLUGIN_init( int ncall )
00106 {
00107    int ii ;
00108    PLUGIN_interface * plint ;     /* will be the output of this routine */
00109 
00110    if( ncall > 3 ) return NULL ;  /* generate interfaces for ncall 0-3 */
00111 
00112    if( ncall == 1 ) return TSGEN_init() ;  /* interface # 1 */
00113    if( ncall == 2 ) return EXP0D_init() ;  /* interface # 2 */
00114 #ifdef ALLOW_LOMO
00115    if( ncall == 3 ) return LOMOR_init() ;  /* interface # 3 */
00116 #else
00117    if( ncall == 3 ) return NULL ;
00118 #endif
00119 
00120    /***** otherwise, do interface # 0 *****/
00121 
00122    /*---------------- set titles and call point ----------------*/
00123 
00124    plint = PLUTO_new_interface( "LSqFit & Dtr" ,
00125                                 "Control LSqFit and LSqDtr Functions" ,
00126                                 helpstring ,
00127                                 PLUGIN_CALL_VIA_MENU , LSQ_main ) ;
00128 
00129    global_plint = plint ;  /* make global copy */
00130 
00131    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00132 
00133    PLUTO_add_hint( plint , "Control LSqFit and LSqDtr Functions" ) ;
00134 
00135    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00136 
00137    /*----- Parameters -----*/
00138 
00139    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00140 
00141    PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
00142 
00143    PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
00144 
00145    /*----- Sinusoid -----*/
00146 
00147    for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
00148       PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00149       PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00150       PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00151    }
00152 
00153    /*----- Timeseries -----*/
00154 
00155    for( ii=0 ; ii < NRMAX_TS ; ii++ ){
00156       PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
00157       PLUTO_add_timeseries( plint , "File" ) ;
00158    }
00159 
00160    /*--------- done with interface setup ---------*/
00161 
00162    PLUTO_register_1D_funcstr( "LSqFit" , LSQ_fitter ) ;
00163    PLUTO_register_1D_funcstr( "LSqDtr" , LSQ_detrend ) ;
00164 
00165    return plint ;
00166 }
00167 
00168 /***************************************************************************
00169   Main routine for this plugin (will be called from AFNI).
00170   If the return string is not NULL, some error transpired, and
00171   AFNI will popup the return string in a message box.
00172 ****************************************************************************/
00173 
00174 char * LSQ_main( PLUGIN_interface * plint )
00175 {
00176    char * str ;
00177    int  ii ;
00178    float * tsar ;
00179 
00180    /*--------- go to first input line ---------*/
00181 
00182    PLUTO_next_option(plint) ;
00183 
00184    str    = PLUTO_get_string(plint) ;
00185    polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
00186 
00187    ignore = PLUTO_get_number(plint) ;
00188 
00189    /*------ loop over remaining options, check their tags, process them -----*/
00190 
00191    nrsin = nrts = 0 ;
00192    do {
00193       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00194 
00195       if( strcmp(str,"Sinusoid") == 0 ){
00196 
00197          sinper[nrsin]  = PLUTO_get_number(plint) ;
00198          sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
00199          if( sinper[nrsin] <= 0.0 )
00200             return "************************\n"
00201                    "Illegal Sinusoid Period!\n"
00202                    "************************"  ;
00203 
00204          nrsin++ ;
00205 
00206       } else if( strcmp(str,"Timeseries") == 0 ){
00207 
00208          tsim[nrts] = PLUTO_get_timeseries(plint) ;
00209 
00210          if( tsim[nrts] == NULL || tsim[nrts]->nx < 3 || tsim[nrts]->kind != MRI_float )
00211             return "*************************\n"
00212                    "Illegal Timeseries Input!\n"
00213                    "*************************"  ;
00214 
00215          tsar = MRI_FLOAT_PTR(tsim[nrts]) ;
00216          for( ii=ignore ; ii < tsim[nrts]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
00217          ignore = ii ;
00218          nrts++ ;
00219 
00220       } else {
00221          return "************************\n"
00222                 "Illegal optiontag found!\n"
00223                 "************************"  ;
00224       }
00225    } while(1) ;
00226 
00227    /*--- nothing left to do until data arrives ---*/
00228 
00229    initialize = 1 ;  /* force re-initialization */
00230 
00231    /*** compute how many ref functions are ordered ***/
00232 
00233    { int nref , ks ;
00234      char str[64] ;
00235 
00236      nref = (polort+1) + nrts ;
00237      for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00238      sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
00239      PLUTO_popup_transient( plint , str ) ;
00240    }
00241 
00242    return NULL ;
00243 }
00244 
00245 /*---------------------------------------------------------------*/
00246 
00247 /** 22 Apr 1997: added label that will go to graphs **/
00248 
00249 void LSQ_fitter( int nt , double to , double dt , float * vec , char ** label )
00250 {
00251    LSQ_worker( nt , dt , vec , TRUE , label ) ;
00252    return ;
00253 }
00254 
00255 void LSQ_detrend( int nt , double to , double dt , float * vec , char ** label )
00256 {
00257    LSQ_worker( nt , dt , vec , FALSE , label ) ;
00258    return ;
00259 }
00260 
00261 static char lbuf[4096] ;  /* 22 Apr 1997: will hold label for graphs */
00262 static char sbuf[256] ;
00263 
00264 void LSQ_worker( int nt , double dt , float * vec , int dofit , char ** label )
00265 {
00266    int nlen , nref ;
00267 
00268    static int nlen_old = -666 , nref_old = -666 ;
00269    static double dt_old = -666.666 ;
00270    static float ** ref = NULL ;
00271    static double * dch = NULL ;
00272 
00273    int ir , ii , ks,jh , j;
00274    float fac , tm , val ;
00275    float * fit , * tsar ;
00276 
00277    /*** compute how many ref functions are ordered ***/
00278 
00279     nref = (polort+1) + nrts ;
00280     for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00281 
00282    /*** do nothing if not enough data to fit ***/
00283 
00284    nlen = nt - ignore ;
00285 
00286    if( nlen <= nref ) return ;  /* do nothing if not enough data to fit */
00287 
00288    /** if data vectors are new length,
00289        or have a new number of reference vectors,
00290        or have a new time step and need sinusoids,
00291        or the initialize flag is set,
00292        then reinitialize reference vectors and Choleski factor **/
00293 
00294    if( nlen != nlen_old || nref != nref_old ||
00295        initialize       || (dt != dt_old && nrsin > 0) ){
00296 
00297       /* free old storage */
00298 
00299       if( ref != NULL ){
00300          for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
00301          free(ref) ;
00302       }
00303       if( dch != NULL ) free(dch) ;
00304 
00305       /* make space for ref vectors */
00306 
00307       ref = (float **) malloc( sizeof(float *) * nref ) ;
00308       if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00309          return;
00310          /* EXIT(1); */
00311       }
00312       for( ir=0 ; ir < nref ; ir++ ){
00313          ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
00314          if( ref[ir] == NULL ){
00315             fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00316             for(j=0;j<=ir;j++)
00317               free(ref[j]);
00318             free(ref);
00319             ref = NULL;
00320             return;
00321             /* EXIT(1); */
00322          }
00323       }
00324       nlen_old = nlen ;
00325       nref_old = nref ;
00326       dt_old   = dt ;
00327 
00328       /**** fill ref vectors ****/
00329 
00330       /* r(t) = 1 */
00331 
00332       for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
00333 
00334       ir = 1 ;
00335       if( polort > 0 ){
00336 
00337          /* r(t) = t - tmid */
00338 
00339          tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
00340          for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
00341          ir = 2 ;
00342 
00343          /* r(t) = (t-tmid)**ir */
00344 
00345          for( ; ir <= polort ; ir++ )
00346             for( ii=0 ; ii < nlen ; ii++ )
00347                ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
00348       }
00349 
00350       if( dt == 0.0 ) dt = 1.0 ;
00351 
00352       /* r(t) = sinusoids */
00353 
00354       for( ks=0 ; ks < nrsin ; ks++ ){
00355          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00356             fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
00357             for( ii=0 ; ii < nlen ; ii++ ){
00358                ref[ir]  [ii] = cos( fac * ii ) ;
00359                ref[ir+1][ii] = sin( fac * ii ) ;
00360             }
00361             ir += 2 ;
00362          }
00363       }
00364 
00365       /* r(t) = timeseries files */
00366 
00367       for( ks=0 ; ks < nrts ; ks++ ){
00368          if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
00369             initialize = 1 ;
00370             fprintf(stderr,"Inadequate time series #%d in LSQ plugin\n\a",ks+1) ;
00371             return ;
00372          }
00373          tsar = MRI_FLOAT_PTR(tsim[ks]) ;
00374          for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = tsar[ii+ignore] ;
00375          ir++ ;
00376       }
00377 
00378       /* Cholesky-ize */
00379 
00380       dch = startup_lsqfit( nlen , NULL , nref , ref ) ;
00381       if( dch == NULL ){
00382          initialize = 1 ;
00383          fprintf(stderr,"Choleski error in LSQ plugin\n\a") ;
00384          return ;
00385       }
00386 
00387       initialize = 0 ;
00388    }
00389 
00390    /** find least squares fit coefficients **/
00391 
00392    fit = delayed_lsqfit( nlen , vec+ignore , nref , ref , dch ) ;
00393 
00394    for( ii=0 ; ii < nlen ; ii++ ){
00395       val = 0.0 ;
00396       for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
00397 
00398       vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
00399    }
00400 
00401    /** 22 Apr 1997: create label if desired by AFNI         **/
00402    /** [This is in static storage, since AFNI will copy it] **/
00403 
00404    if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00405 
00406       lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00407 
00408       /** for each reference, make a string into sbuf **/
00409 
00410       ir = 0 ;
00411       sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
00412       strcat(lbuf,sbuf) ;
00413 
00414       for( ; ir <= polort ; ){
00415          sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir++] ) ;
00416          strcat(lbuf,sbuf) ;
00417       }
00418 
00419       for( ks=0 ; ks < nrsin ; ks++ ){
00420          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00421             fac = sinper[ks] / (jh+1) ;
00422             sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00423             strcat(lbuf,sbuf) ;
00424             sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00425             strcat(lbuf,sbuf) ;
00426          }
00427       }
00428 
00429       for( ks=0 ; ks < nrts ; ks++ ){
00430          sprintf(sbuf,"Coef of %s = %g\n" , tsim[ks]->name , fit[ir++] ) ;
00431          strcat(lbuf,sbuf) ;
00432       }
00433 
00434       *label = lbuf ;  /* send address of lbuf back in what label points to */
00435    }
00436 
00437    free(fit) ;
00438    return ;
00439 }
00440 
00441 /*********************************************************************************/
00442 
00443 PLUGIN_interface * TSGEN_init(void)
00444 {
00445    PLUGIN_interface * plint ;
00446 
00447    /*---------------- set titles and call point ----------------*/
00448 
00449    plint = PLUTO_new_interface( "TS Generate" ,
00450                                 "Generate a Timeseries" ,
00451                                 plehstring ,
00452                                 PLUGIN_CALL_VIA_MENU , TSGEN_main ) ;
00453 
00454    PLUTO_add_hint( plint , "Generate a 1D Timeseries" ) ;
00455 
00456    /*----- Parameters -----*/
00457 
00458    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00459    PLUTO_add_number( plint , "Delta"  , 0,99999,1, 0 , TRUE ) ;
00460    PLUTO_add_number( plint , "Length" , 3,9999,0,3   , TRUE ) ;
00461 
00462    /*----- Output -----*/
00463 
00464    PLUTO_add_option( plint , "Output" , "Output" , TRUE ) ;
00465    PLUTO_add_string( plint , "Label" , 0,NULL , 19 ) ;
00466 
00467    /*----- Polynomial -----*/
00468 
00469    PLUTO_add_option( plint , "Polynomial" , "Polynomial" , FALSE ) ;
00470    PLUTO_add_number( plint , "Order" , 2,20,0,2 , FALSE ) ;
00471 
00472    /*----- Sinusoid -----*/
00473 
00474    PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00475    PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00476    PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00477 
00478    return plint ;
00479 }
00480 
00481 char * TSGEN_main( PLUGIN_interface * plint )
00482 {
00483    char * label , * str ;
00484    int  ii , jj ;
00485    float * tsar ;
00486    MRI_IMAGE * tsim ;
00487    float delta , period=0.0 ;
00488    int   nx , ny=0 , npol=0 , nharm=-1 ;
00489    int pp ;
00490    double fac , val ;
00491 
00492    /*--------- go to first input line ---------*/
00493 
00494    PLUTO_next_option(plint) ;
00495 
00496    delta = PLUTO_get_number(plint) ;
00497    if( delta <= 0.0 ) return "**********************\n"
00498                              "Illegal value of Delta\n"
00499                              "**********************"  ;
00500 
00501    nx = PLUTO_get_number(plint) ;
00502 
00503    /*----- next input line -----*/
00504 
00505    PLUTO_next_option(plint) ;
00506    label = PLUTO_get_string(plint) ;
00507    if( label == NULL || strlen(label) == 0 ) return "**********************\n"
00508                                                     "Illegal value of Label\n"
00509                                                     "**********************"  ;
00510 
00511    /*----- rest of input lines -----*/
00512 
00513    do {
00514       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00515 
00516       if( strcmp(str,"Sinusoid") == 0 ){
00517 
00518          period = PLUTO_get_number(plint) ;
00519          nharm  = PLUTO_get_number(plint) - 1.0 ;
00520 
00521          if( period <= 0.0 ) return "***********************\n"
00522                                     "Illegal Sinusoid Period\n"
00523                                     "***********************"  ;
00524 
00525 
00526       } else if( strcmp(str,"Polynomial") == 0 ){
00527 
00528          npol = PLUTO_get_number(plint) ;
00529 
00530       } else {
00531          return "***********************\n"
00532                 "Illegal optiontag found\n"
00533                 "***********************"  ;
00534       }
00535    } while(1) ;
00536 
00537    /********** Make the timeseries ***********/
00538 
00539    ny = 0 ;
00540    if( npol > 0 )     ny  = npol-1 ;
00541    if( period > 0.0 ) ny += 2*(nharm+1) ;
00542 
00543    if( ny < 1 ) return "***********************\n"
00544                        "No timeseries specified\n"
00545                        "***********************"  ;
00546 
00547    tsim = mri_new( nx , ny , MRI_float ) ;
00548    jj   = 0 ;
00549 
00550    fac  = 1.99999 / (nx-1) ;
00551    for( pp=2 ; pp <= npol ; pp++,jj++ ){
00552 
00553       tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
00554 
00555       for( ii=0 ; ii < nx ; ii++ ){
00556          val = fac * ii - 0.999995 ;
00557          tsar[ii] = cos( pp * acos(val) ) ;
00558       }
00559    }
00560 
00561    for( pp=0 ; pp <= nharm ; pp++ , jj+=2 ){
00562       fac  = (2.0*PI) * delta * (pp+1) / period ;
00563       tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
00564 
00565       for( ii=0 ; ii < nx ; ii++ ){
00566          tsar[ii]    = cos( fac * ii ) ;
00567          tsar[ii+nx] = sin( fac * ii ) ;
00568       }
00569    }
00570 
00571    PLUTO_register_timeseries( label , tsim ) ;
00572    mri_free(tsim) ;
00573    return NULL ;
00574 }
00575 
00576 /***************************************************************************/
00577 
00578 #include "parser.h"
00579 
00580 static char fredstring[] =
00581   " Purpose: Control the Expr 0D Transformation function\n"
00582   "\n"
00583   " Variable   = letter used as input to expression\n"
00584   " Expression = arithmetic expression to evaluate\n"
00585 ;
00586 
00587 #define NALPHA 26
00588 static char * vstring[NALPHA] = {
00589   " A ",  " B ",  " C ",  " D ",  " E ",
00590   " F ",  " G ",  " H ",  " I ",  " J ",
00591   " K ",  " L ",  " M ",  " N ",  " O ",
00592   " P ",  " Q ",  " R ",  " S ",  " T ",
00593   " U ",  " V ",  " W ",  " X ",  " Y ",  " Z " } ;
00594 
00595 static int exp0d_var = 23 ;
00596 
00597 static PARSER_code * exp0d_pc = NULL ;
00598 
00599 static PLUGIN_interface *plint_EXP0D=NULL ;
00600 
00601 static void EXP0D_func_init(void)   /* 21 Jul 2003 */
00602 {
00603    PLUG_startup_plugin_CB( NULL , (XtPointer)plint_EXP0D , NULL ) ;
00604 }
00605 
00606 
00607 PLUGIN_interface * EXP0D_init(void)
00608 {
00609    PLUGIN_interface * plint ;
00610 
00611    /*---------------- set titles and call point ----------------*/
00612 
00613    plint = PLUTO_new_interface( "Expr 0D" ,
00614                                 "Control the Expr 0D transformation" ,
00615                                 fredstring ,
00616                                 PLUGIN_CALL_VIA_MENU , EXP0D_main ) ;
00617 
00618    PLUTO_add_option( plint , "Variable" , "Variable" , TRUE ) ;
00619    PLUTO_add_string( plint , NULL , NALPHA,vstring , exp0d_var ) ;
00620 
00621    PLUTO_add_option( plint , "Expression" , "Expression" , TRUE ) ;
00622    PLUTO_add_string( plint , NULL , 0,NULL , 50 ) ;
00623 
00624    PLUTO_register_0D_function( "Expr 0D" , EXP0D_worker ) ;
00625 
00626    plint_EXP0D = plint ;
00627    AFNI_register_nD_func_init( 0 , (generic_func *)EXP0D_func_init ) ;  /* 21 Jul 2003 */
00628 
00629    return plint ;
00630 }
00631 
00632 char * EXP0D_main( PLUGIN_interface * plint )
00633 {
00634    char * str ;
00635 
00636    PLUTO_next_option(plint) ;
00637    str       = PLUTO_get_string(plint) ;
00638    exp0d_var = PLUTO_string_index( str , NALPHA , vstring ) ;
00639 
00640    if( exp0d_pc != NULL ){ free(exp0d_pc) ; exp0d_pc = NULL ; }
00641    PLUTO_next_option(plint) ;
00642    str       = PLUTO_get_string(plint) ;
00643    exp0d_pc  = PARSER_generate_code( str ) ;
00644 
00645    if( exp0d_pc == NULL ) return "*******************************\n"
00646                                  "Error when compiling expression\n"
00647                                  "*******************************"  ;
00648 
00649    return NULL ;
00650 }
00651 
00652 #define VSIZE 64
00653 
00654 void EXP0D_worker( int num , float * vec )
00655 {
00656    int ii , jj , jbot,jtop ;
00657 
00658    static int first = 1 ;
00659    static double * atoz[26] ;
00660    static double   tvec[VSIZE] ;
00661 
00662    if( num <= 0 || vec == NULL || exp0d_pc == NULL ) return ;
00663 
00664 #if 0
00665 fprintf(stderr,"Enter EXP0D_worker\n") ;
00666 #endif
00667 
00668    if( first ){
00669       for( ii=0 ; ii < 26 ; ii++)
00670         atoz[ii] = (double *) malloc(sizeof(double) * VSIZE ) ;
00671       first = 0 ;
00672 #if 0
00673 fprintf(stderr,"Allocated atoz\n") ;
00674 #endif
00675    }
00676 
00677    for( ii=0 ; ii < 26 ; ii++ )
00678       for (jj=0; jj<VSIZE; jj++) atoz[ii][jj] = 0.0 ;
00679 
00680 #if 0
00681 fprintf(stderr,"Zeroed atoz\n") ;
00682 #endif
00683 
00684    for( ii=0 ; ii < num ; ii+=VSIZE ){
00685       jbot = ii ;
00686       jtop = MIN( ii + VSIZE , num ) ;
00687 
00688       for( jj=jbot ; jj < jtop ; jj ++ ) atoz[exp0d_var][jj-ii] = vec[jj] ;
00689 
00690       PARSER_evaluate_vector( exp0d_pc , atoz , jtop-jbot , tvec ) ;
00691 
00692       for( jj=jbot ; jj < jtop ; jj ++ ) vec[jj] = tvec[jj-ii] ;
00693    }
00694 
00695 #if 0
00696 fprintf(stderr,"Exit EXP0D_worker\n") ;
00697 #endif
00698 
00699    return ;
00700 }
00701 
00702 /*------------------------------------------------------------------*/
00703 
00704 #ifdef ALLOW_LOMO
00705 static int lomo_order  = 3 ;
00706 static int lomo_levels = 32 ;
00707 int lomo_regress( int , int , int * , int * ) ;
00708 
00709 PLUGIN_interface * LOMOR_init(void)
00710 {
00711    PLUGIN_interface * plint ;
00712 
00713    /*---------------- set titles and call point ----------------*/
00714 
00715    plint = PLUTO_new_interface( "Lomo Regression" ,
00716                                 "Control the Local Monotone Regression" ,
00717                                 NULL ,
00718                                 PLUGIN_CALL_VIA_MENU , LOMOR_main ) ;
00719 
00720    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00721    PLUTO_add_number( plint , "Order" , 3, 20,0,lomo_order  , FALSE ) ;
00722    PLUTO_add_number( plint , "Levels", 4,128,0,lomo_levels , FALSE ) ;
00723 
00724    PLUTO_register_1D_function( "Lomo 1D" , LOMOR_fitter ) ;
00725 
00726    return plint ;
00727 }
00728 
00729 char * LOMOR_main( PLUGIN_interface * plint )
00730 {
00731    PLUTO_next_option(plint) ;
00732    lomo_order  = PLUTO_get_number(plint) ;
00733    lomo_levels = PLUTO_get_number(plint) ;
00734    return NULL ;
00735 }
00736 
00737 static int lnum   = 0 ;
00738 static int * xin  = NULL ;
00739 static int * yout = NULL ;
00740 
00741 void LOMOR_fitter( int num , double to , double dt , float * vec )
00742 {
00743    int ii ;
00744    float top , bot , scl ;
00745 
00746    if( num <= 3 || vec == NULL ) return ;
00747 
00748    if( lnum < num ){
00749       if( xin != NULL ){ free(xin) ; free(yout) ; }
00750       xin  = (int *) malloc( sizeof(int) * num ) ;
00751       yout = (int *) malloc( sizeof(int) * num ) ;
00752       if( xin == NULL || yout == NULL ) return ;
00753       lnum = num ;
00754    }
00755 
00756    bot = top = vec[0] ;
00757    for( ii=1 ; ii < num ; ii++ ){
00758            if( vec[ii] < bot ) bot = vec[ii] ;
00759       else if( vec[ii] > top ) top = vec[ii] ;
00760    }
00761 
00762    if( bot >= top ) return ;
00763    scl = (lomo_levels - 0.01) / (top-bot) ;
00764 
00765 fprintf(stderr,"LOMO: range %f .. %f; scl=%f\n",bot,top,scl) ;
00766 
00767    for( ii=0 ; ii < num ; ii++ ) xin[ii] = (int)((vec[ii]-bot)*scl) ;
00768 
00769    ii = lomo_regress( num , lomo_order , xin , yout ) ;
00770    if( ii == -1 ) return ;
00771 
00772    scl = 1.0 / scl ;
00773    for( ii=0 ; ii < num ; ii++ ) vec[ii] = bot + yout[ii]*scl ;
00774    return ;
00775 }
00776 
00777 /*******************************************************/
00778 /* Fast Digital Locally Monotonic Regression           */
00779 /*******************************************************/
00780 /*               Copyright (c) 1995                    */
00781 /*    University of Maryland at College Park           */
00782 /*               All Rights Reserved                   */
00783 /*******************************************************/
00784 /*  by Nicholas Sidiropoulos, Aug.  1995               */
00785 /*******************************************************/
00786 /* compile using -lm option                            */
00787 /*******************************************************/
00788 /* input: from stdinput.dat: standard matlab vector    */
00789 /*        (i.e., ASCII file containing                 */
00790 /*               a long line of N vector elements      */
00791 /*               separated by spaces)                  */
00792 /* output: in stdoutput.dat.*: same format as input    */
00793 /* control: in stdcontrol_switch.dat: the              */
00794 /*          ``effective'' M is  */
00795 /*          the first entry here, followed by `\n`     */
00796 /*          Rest: fixed to 1 (future option)           */
00797 /*          So stdcontrol_switch .dat should be ASCII  */
00798 /*          file starting with the effective M \n      */
00799 /*          and followed by N `1''s \n, e.g., for M=10 */
00800 /*  10<newline>1<newline>...1<newline>                 */
00801 /*                 a total of N ones                   */
00802 /*******************************************************/
00803 /* This is just an ``Academic'' piece of software -    */
00804 /* it has been checked for correctness to the best     */
00805 /* of my ability, however, no guarantees whatsoever    */
00806 /* are given - all disclaimers here!                   */
00807 /* It has been optimized for minimum development effort*/
00808 /* and NOT for optimum speed and/or minimum comput.    */
00809 /* resources. Note, in particular, that I do not take  */
00810 /* advantage of trellis path merging to reduce storage */
00811 /* requirements. As a result, depending on your choice */
00812 /* of parameters M,N,R, the program may require        */
00813 /* considerable amounts of static memory. Therefore,   */
00814 /* if your computer lacks it you need to rlogin to     */
00815 /* a machine (like oxygen or apollo at SIL lab) which  */
00816 /* has sufficient memory                               */
00817 /*******************************************************/
00818 
00819 /*=====================================================*/
00820 /*== Modified Feb 1997 by RWCox, to be a C function. ==*/
00821 /*=====================================================*/
00822 
00823 #include <stdio.h>
00824 #include <math.h>
00825 #include <string.h>
00826 
00827 typedef struct _state {
00828   int value, length, cumcost;
00829   struct _state *prevstate;
00830 } state;
00831 
00832 #ifdef ABS
00833 #undef ABS
00834 #endif
00835 #define ABS(x) (((x)<0) ? (-x) : (x))
00836 
00837 #define USE_STATIC_TRELLIS
00838 #ifdef USE_STATIC_TRELLIS
00839 #  define MMAX 35    /* (strictly) > max lomo degree   */
00840 #  define NMAX 512   /* input data length              */
00841 #  define RMAX 100   /* range of input: 0 -> R-1 inc.  */
00842    static state TRELLIS[RMAX][2*MMAX][NMAX] ;
00843 #  define trellis(i,j,k) TRELLIS[i][j][k]
00844 #else
00845    static state * TRELLIS = NULL ;
00846    static int     TDIM1 , TDIM2 , TDIM12 , TDIM3 ;
00847 #  define trellis(i,j,k) TRELLIS[i+j*TDIM1+k*TDIM12]
00848 #endif
00849 
00850 /*=====================================================*/
00851 /*==   N   = number of input points                  ==*/
00852 /*==   m   = local monotonicity order to impose      ==*/
00853 /*==   yin = pointer to inputs  (array of length N)  ==*/
00854 /*==   x   = pointer to outputs (array of length N)  ==*/
00855 /*==                                                 ==*/
00856 /*==   return value is 0 if all is OK, -1 if not     ==*/
00857 /*=====================================================*/
00858 
00859 int lomo_regress( int N , int m , int * yin , int * x )
00860 {
00861   int base , R , * y ;
00862   int n,v,l,pv,pl,i,j;
00863   int cost, maxcost;
00864   int peakval;
00865   state dummy_state;      /* dummy initial state */
00866 
00867   /*== check inputs for OK-ness ==*/
00868 
00869   if( N < 3 || m < 3 || m >= N || yin == NULL || x == NULL ) return -1 ;
00870 
00871   /*== Compute range of input into R ==*/
00872 
00873   y = (int *) malloc(sizeof(int)*N) ; if( y == NULL ) return -1 ;
00874   base = yin[0] ;
00875   for( i=1 ; i < N ; i++ ) if( yin[i] < base ) base = yin[i] ;
00876   R = y[0] = yin[0] - base ;
00877   for( i=1 ; i < N ; i++ ){
00878      y[i] = yin[i] - base ; if( y[i] > R ) R = y[i] ;
00879   }
00880   R++ ;
00881   if( R == 1 ){
00882       for( i=0 ; i < N ; i++ ) x[i] = yin[i] ;
00883       free(y) ; return 0 ;
00884   }
00885 
00886 fprintf(stderr,"LOMO: %d points; %d levels; %d order\n",N,R,m) ;
00887 
00888   /* compute maxcost */
00889 
00890   maxcost = 0;
00891   for (n = 0; n < N; n++) maxcost += ABS(y[n]);
00892 
00893   /* now init FLOMOR : */
00894 
00895   dummy_state.value   = -1; dummy_state.length    = m   ;
00896   dummy_state.cumcost = 0 ; dummy_state.prevstate = NULL;
00897 
00898 #ifndef USE_STATIC_TRELLIS
00899   /*== malloc trellis space ==*/
00900 
00901   TDIM1 = R ; TDIM2 = 2*m+2 ; TDIM12 = TDIM1*TDIM2 ; TDIM3 = N ;
00902   TRELLIS = (state *) calloc( TDIM1*TDIM2*TDIM3 , sizeof(state) ) ;
00903   if( TRELLIS == NULL ){ free(y) ; return -1 ; }
00904 #endif
00905 
00906 fprintf(stderr,"LOMO: init trellis\n") ;
00907 
00908   for (n = 0; n < N; n++)
00909   {
00910     for (l = 1; l <= 2*m; l++)
00911     {
00912       for (v = 0; v < R; v++) { trellis(v,l,n).value = v;
00913                                 trellis(v,l,n).length = l;
00914                                 if ((n == 0) && ((l == 1) || (l == m+1)))
00915                                 { trellis(v,l,n).cumcost   = ABS((v-y[0]));
00916                                   trellis(v,l,n).prevstate = &dummy_state;
00917                                 }
00918                                 else
00919                                 { trellis(v,l,n).cumcost   = maxcost;
00920                                   trellis(v,l,n).prevstate = NULL;
00921                                 }
00922                               }
00923     }
00924   }
00925 
00926 /************************ main FLOMOR loop ************************************/
00927 /* note that l = -1, ..., -m is mapped to l = m+1, ..., 2m resp.              */
00928 /******************************************************************************/
00929 
00930 fprintf(stderr,"LOMO: compute trellis") ; fflush(stderr) ;
00931 
00932   for (n = 1; n < N; n++)
00933   {
00934     /* states of the first type: (v,1) */
00935 
00936       fprintf(stderr,".") ; fflush(stderr) ;
00937 
00938       for (v = 0; v < R; v++)
00939       {
00940         for (pv = 0; pv < v; pv++)
00941         {
00942           for (pl = 1; pl <= m; pl++)
00943           {
00944             if (trellis(pv,pl,n-1).cumcost != maxcost)
00945             {
00946               cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
00947               if (cost < trellis(v,1,n).cumcost)
00948               {
00949                 trellis(v,1,n).cumcost = cost;
00950                 trellis(v,1,n).prevstate = &trellis(pv,pl,n-1);
00951               }
00952             }
00953           }
00954           if (trellis(pv,2*m,n-1).cumcost != maxcost)
00955           {
00956             cost = trellis(pv,2*m,n-1).cumcost + ABS((v-y[n]));
00957             if (cost < trellis(v,1,n).cumcost)
00958             {
00959               trellis(v,1,n).cumcost = cost;
00960               trellis(v,1,n).prevstate = &trellis(pv,2*m,n-1);
00961             }
00962           }
00963         }
00964       }
00965 
00966     /* states of the second type: (v,-1) */
00967 
00968       for (v = 0; v < R; v++)
00969       {
00970         for (pv = R - 1; pv > v; pv--)
00971         {
00972           for (pl = m + 1; pl <= 2*m; pl++)
00973           {
00974             if (trellis(pv,pl,n-1).cumcost != maxcost)
00975             {
00976               cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
00977               if (cost < trellis(v,m+1,n).cumcost)
00978               {
00979                 trellis(v,m+1,n).cumcost = cost;
00980                 trellis(v,m+1,n).prevstate = &trellis(pv,pl,n-1);
00981               }
00982             }
00983           }
00984           if (trellis(pv,m,n-1).cumcost != maxcost)
00985           {
00986             cost = trellis(pv,m,n-1).cumcost + ABS((v-y[n]));
00987             if (cost < trellis(v,m+1,n).cumcost)
00988             {
00989               trellis(v,m+1,n).cumcost = cost;
00990               trellis(v,m+1,n).prevstate = &trellis(pv,m,n-1);
00991             }
00992           }
00993         }
00994       }
00995 
00996     /* states of the third type: (v,l), 1 < l < m */
00997 
00998     for (l = 2; l < m; l++)
00999     {
01000       for (v = 0; v < R; v++)
01001       {
01002         if (trellis(v,l-1,n-1).cumcost != maxcost)
01003         {
01004           trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
01005           trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
01006         }
01007       }
01008     }
01009 
01010     /* states of the fourth type: (v,l), -m < l < -1 */
01011 
01012     for (l = m+2; l < 2*m; l++)
01013     {
01014       for (v = 0; v < R; v++)
01015       {
01016         if (trellis(v,l-1,n-1).cumcost != maxcost)
01017         {
01018           trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
01019           trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
01020         }
01021       }
01022     }
01023 
01024     /* states of the fifth type: (v,l), l = m */
01025 
01026     for (v = 0; v < R; v++)
01027     {
01028       if (trellis(v,m,n-1).cumcost != maxcost)
01029       {
01030         cost = trellis(v,m,n-1).cumcost + ABS((v-y[n]));
01031         if (cost < trellis(v,m,n).cumcost)
01032         {
01033           trellis(v,m,n).cumcost = cost;
01034           trellis(v,m,n).prevstate = &trellis(v,m,n-1);
01035         }
01036       }
01037 
01038       if (trellis(v,m-1,n-1).cumcost != maxcost)
01039       {
01040         cost = trellis(v,m-1,n-1).cumcost + ABS((v-y[n]));
01041         if (cost < trellis(v,m,n).cumcost)
01042         {
01043           trellis(v,m,n).cumcost = cost;
01044           trellis(v,m,n).prevstate = &trellis(v,m-1,n-1);
01045         }
01046       }
01047     }
01048 
01049     /* states of the sixth type: (v,l), l = -m */
01050 
01051     for (v = 0; v < R; v++)
01052     {
01053       if (trellis(v,2*m,n-1).cumcost != maxcost)
01054       {
01055         cost = trellis(v,2*m,n-1).cumcost + ABS((v-y[n]));
01056         if (cost < trellis(v,2*m,n).cumcost)
01057         {
01058           trellis(v,2*m,n).cumcost = cost;
01059           trellis(v,2*m,n).prevstate = &trellis(v,2*m,n-1);
01060         }
01061       }
01062 
01063       if (trellis(v,2*m-1,n-1).cumcost != maxcost)
01064       {
01065         cost = trellis(v,2*m-1,n-1).cumcost + ABS((v-y[n]));
01066         if (cost < trellis(v,2*m,n).cumcost)
01067         {
01068           trellis(v,2*m,n).cumcost = cost;
01069           trellis(v,2*m,n).prevstate = &trellis(v,2*m-1,n-1);
01070         }
01071       }
01072     }
01073 
01074 
01075   }
01076 
01077 /************************ eof main FLOMOR loop *******************************/
01078 
01079   /* now pick best path, and write it to x[n] */
01080 
01081 fprintf(stderr,"\nLOMO: pick best path\n") ; fflush(stderr) ;
01082 
01083   cost = maxcost;
01084   for (l = 1; l <= m; l++)
01085   {
01086     for (v = 0; v < R; v++)
01087     {
01088       if (trellis(v,l,N-1).cumcost < cost)
01089       {
01090         cost = trellis(v,l,N-1).cumcost;
01091         pl = l; pv = v;
01092       }
01093     }
01094   }
01095 
01096   /* now traverse backwards */
01097 
01098   if (cost < maxcost)
01099   {
01100     x[N-1] = pv;
01101 fprintf(stderr,"  [%d] = %d\n",N-1,pv) ;
01102 
01103     for (n = N-2; n >=0; n--)
01104     {
01105 fprintf(stderr,"  [%d] = %d",n,trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
01106       x[n] = trellis(pv,pl,n+1).prevstate->value;
01107 fprintf(stderr,"  pv = %d",trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
01108       pv = trellis(pv,pl,n+1).prevstate->value;
01109 fprintf(stderr,"  pl = %d",trellis(pv,pl,n+1).prevstate->length) ; fflush(stderr) ;
01110       pl = trellis(pv,pl,n+1).prevstate->length;
01111 fprintf(stderr,"\n") ;
01112     }
01113   }
01114   else { for (n=0;n<N;n++) x[n] = 0; /* as good as any... */ }
01115 
01116   /*== add base back to output ==*/
01117 
01118 fprintf(stderr,"\nLOMO: add base back on\n") ;
01119 
01120   for( i=0 ; i < N ; i++ ) x[n] += base ;
01121 
01122 fprintf(stderr,"LOMO: exit regression\n") ;
01123 
01124   free(y) ;
01125 #ifndef USE_STATIC_TRELLIS
01126   free(TRELLIS) ;
01127 #endif
01128 
01129   return 0 ;
01130 }
01131 #endif /* ALLOW_LOMO */