00001 #include "afni.h"
00002 
00003 #ifndef ALLOW_PLUGINS
00004 #  error "Plugins not properly set up -- see machdep.h"
00005 #endif
00006 
00007 /***********************************************************************
00008   Plugin to provide a L1 fitting 1D function for graphs
00009 ************************************************************************/
00010 
00011 /*------------- string to 'help' the user -------------*/
00012 
00013 static char helpstring[] =
00014    " Purpose: Control the 'L1_Fit' and 'L1_Dtr 1D functions.\n"
00015    "\n"
00016    " Parameters:  Baseline = 'Constant' or 'Linear'\n"
00017    "                           Is the baseline 'a' or 'a+b*t'?\n"
00018    "              Ignore   = Number of points to ignore at\n"
00019    "                           start of each timeseries.\n"
00020    " \n"
00021    " Sinusoids:   Period    = Fundamental period to use.\n"
00022    "              Harmonics = Number of overtones to use.\n"
00023    " \n"
00024    " Timeseries:  File      = Input timeseries file to use.\n"
00025 ;
00026 
00027 /*------- Strings for baseline control ------*/
00028 
00029 #define NBASE 4
00030 static char * baseline_strings[NBASE] = { "Constant", "Linear", "Quadratic", "Cubic" } ;
00031 
00032 /*--------------- prototypes for internal routines ---------------*/
00033 
00034 char * L1F_main( PLUGIN_interface * ) ;  /* the entry point */
00035 
00036 void L1F_fitter ( int nt, double to, double dt, float * vec, char ** label ) ;
00037 void L1F_detrend( int nt, double to, double dt, float * vec, char ** label ) ;
00038 void L1F_worker ( int nt, double dt, float * vec, int dofit, char ** label ) ;
00039 
00040 /*---------------- global data -------------------*/
00041 
00042 static PLUGIN_interface * global_plint = NULL ;
00043 
00044 #define NRMAX_SIN 2
00045 #define NRMAX_TS  2
00046 #define HARM_MAX  22
00047 
00048 static int polort=1 , ignore=3 , nrsin=0 , nrts=0 , initialize=1 ;
00049 static float sinper[NRMAX_SIN] ;
00050 static int   sinharm[NRMAX_SIN] ;
00051 static MRI_IMAGE * tsim[NRMAX_TS] ;
00052 
00053 /***********************************************************************
00054    Set up the interface to the user:
00055     1) Create a new interface using "PLUTO_new_interface";
00056 
00057     2) For each line of inputs, create the line with "PLUTO_add_option"
00058          (this line of inputs can be optional or mandatory);
00059 
00060     3) For each item on the line, create the item with
00061         "PLUTO_add_dataset"    for a dataset chooser,
00062         "PLUTO_add_string"     for a string chooser,
00063         "PLUTO_add_number"     for a number chooser,
00064         "PLUTO_add_timeseries" for a timeseries chooser.
00065 ************************************************************************/
00066 
00067 
00068 DEFINE_PLUGIN_PROTOTYPE
00069 
00070 PLUGIN_interface * PLUGIN_init( int ncall )
00071 {
00072    int ii ;
00073    PLUGIN_interface * plint ;     /* will be the output of this routine */
00074 
00075    if( ncall > 0 ) return NULL ;
00076 
00077    /*---------------- set titles and call point ----------------*/
00078 
00079    plint = PLUTO_new_interface( "L1_Fit & Dtr" ,
00080                                 "Control L1_Fit and L1_Dtr Functions" ,
00081                                 helpstring ,
00082                                 PLUGIN_CALL_VIA_MENU , L1F_main ) ;
00083 
00084    global_plint = plint ;  /* make global copy */
00085 
00086    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00087 
00088    PLUTO_add_hint( plint , "Control L1_Fit and L1_Dtr Functions" ) ;
00089 
00090    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00091 
00092    /*----- Parameters -----*/
00093 
00094    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00095 
00096    PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
00097 
00098    PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
00099 
00100    /*----- Sinusoid -----*/
00101 
00102    for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
00103       PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00104       PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00105       PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00106    }
00107 
00108    /*----- Timeseries -----*/
00109 
00110    for( ii=0 ; ii < NRMAX_TS ; ii++ ){
00111       PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
00112       PLUTO_add_timeseries( plint , "File" ) ;
00113    }
00114 
00115    /*--------- done with interface setup ---------*/
00116 
00117    PLUTO_register_1D_funcstr( "L1_Fit" , L1F_fitter ) ;
00118    PLUTO_register_1D_funcstr( "L1_Dtr" , L1F_detrend ) ;
00119 
00120    return plint ;
00121 }
00122 
00123 /***************************************************************************
00124   Main routine for this plugin (will be called from AFNI).
00125   If the return string is not NULL, some error transpired, and
00126   AFNI will popup the return string in a message box.
00127 ****************************************************************************/
00128 
00129 char * L1F_main( PLUGIN_interface * plint )
00130 {
00131    char * str ;
00132    int  ii ;
00133    float * tsar ;
00134 
00135    /*--------- go to first input line ---------*/
00136 
00137    PLUTO_next_option(plint) ;
00138 
00139    str    = PLUTO_get_string(plint) ;
00140    polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
00141 
00142    ignore = PLUTO_get_number(plint) ;
00143 
00144    /*------ loop over remaining options, check their tags, process them -----*/
00145 
00146    nrsin = nrts = 0 ;
00147    do {
00148       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00149 
00150       if( strcmp(str,"Sinusoid") == 0 ){
00151 
00152          sinper[nrsin]  = PLUTO_get_number(plint) ;
00153          sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
00154          if( sinper[nrsin] <= 0.0 )
00155             return "************************\n"
00156                    "Illegal Sinusoid Period!\n"
00157                    "************************"  ;
00158 
00159          nrsin++ ;
00160 
00161       } else if( strcmp(str,"Timeseries") == 0 ){
00162 
00163          tsim[nrts] = PLUTO_get_timeseries(plint) ;
00164 
00165          if( tsim[nrts] == NULL || tsim[nrts]->nx < 3 || tsim[nrts]->kind != MRI_float )
00166             return "*************************\n"
00167                    "Illegal Timeseries Input!\n"
00168                    "*************************"  ;
00169 
00170          tsar = MRI_FLOAT_PTR(tsim[nrts]) ;
00171          for( ii=ignore ; ii < tsim[nrts]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
00172          ignore = ii ;
00173          nrts++ ;
00174 
00175       } else {
00176          return "************************\n"
00177                 "Illegal optiontag found!\n"
00178                 "************************"  ;
00179       }
00180    } while(1) ;
00181 
00182    /*--- nothing left to do until data arrives ---*/
00183 
00184    initialize = 1 ;  /* force re-initialization */
00185 
00186    /*** compute how many ref functions are ordered ***/
00187 
00188    { int nref , ks ;
00189      char str[64] ;
00190 
00191      nref = (polort+1) + nrts ;
00192      for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00193      sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
00194      PLUTO_popup_transient( plint , str ) ;
00195    }
00196 
00197    return NULL ;
00198 }
00199 
00200 /*---------------------------------------------------------------*/
00201 
00202 /** 22 Apr 1997: added label that will go to graphs **/
00203 
00204 void L1F_fitter( int nt , double to , double dt , float * vec , char ** label )
00205 {
00206    L1F_worker( nt , dt , vec , TRUE , label ) ;
00207    return ;
00208 }
00209 
00210 void L1F_detrend( int nt , double to , double dt , float * vec , char ** label )
00211 {
00212    L1F_worker( nt , dt , vec , FALSE , label ) ;
00213    return ;
00214 }
00215 
00216 static char lbuf[4096] ;  /* 22 Apr 1997: will hold label for graphs */
00217 static char sbuf[256] ;
00218 
00219 void L1F_worker( int nt , double dt , float * vec , int dofit , char ** label )
00220 {
00221    int nlen , nref ;
00222 
00223    static int nlen_old = -666 , nref_old = -666 ;
00224    static double dt_old = -666.666 ;
00225    static float ** ref = NULL ;
00226    static float *  fit = NULL ;
00227 
00228    int ir , ii , ks,jh ;
00229    float fac , tm , val , cls ;
00230    float * tsar ;
00231 
00232    /*** compute how many ref functions are ordered ***/
00233 
00234    nref = (polort+1) + nrts ;
00235    for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00236 
00237    /*** do nothing if not enough data to fit ***/
00238 
00239    nlen = nt - ignore ;
00240 
00241    if( nlen <= nref ) return ;  /* do nothing if not enough data to fit */
00242 
00243    /** if data vectors are new length,
00244        or have a new number of reference vectors,
00245        or have a new time step and need sinusoids,
00246        or the initialize flag is set,
00247        then reinitialize reference vectors and Choleski factor **/
00248 
00249    if( nlen != nlen_old || nref != nref_old ||
00250        initialize       || (dt != dt_old && nrsin > 0) ){
00251 
00252       /* free old storage */
00253 
00254       if( ref != NULL ){
00255          for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
00256          free(ref) ;
00257       }
00258       if( fit != NULL ) free(fit) ;
00259 
00260       /* make space for ref vectors */
00261 
00262       ref = (float **) malloc( sizeof(float *) * nref ) ;
00263       if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00264          return;
00265          /*EXIT(1);*/
00266       }
00267       for( ir=0 ; ir < nref ; ir++ ){
00268          ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
00269          if( ref[ir] == NULL )
00270             {fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00271             return;
00272             /* EXIT(1); */
00273          }
00274       }
00275       nlen_old = nlen ;
00276       nref_old = nref ;
00277       dt_old   = dt ;
00278 
00279       /**** fill ref vectors ****/
00280 
00281       /* r(t) = 1 */
00282 
00283       for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
00284 
00285       ir = 1 ;
00286       if( polort > 0 ){
00287 
00288          /* r(t) = t - tmid */
00289 
00290          tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
00291          for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
00292          ir = 2 ;
00293 
00294          /* r(t) = (t-tmid)**ir */
00295 
00296          for( ; ir <= polort ; ir++ )
00297             for( ii=0 ; ii < nlen ; ii++ )
00298                ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
00299       }
00300 
00301       if( dt == 0.0 ) dt = 1.0 ;
00302 
00303       /* r(t) = sinusoids */
00304 
00305       for( ks=0 ; ks < nrsin ; ks++ ){
00306          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00307             fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
00308             for( ii=0 ; ii < nlen ; ii++ ){
00309                ref[ir]  [ii] = cos( fac * ii ) ;
00310                ref[ir+1][ii] = sin( fac * ii ) ;
00311             }
00312             ir += 2 ;
00313          }
00314       }
00315 
00316       /* r(t) = timeseries files */
00317 
00318       for( ks=0 ; ks < nrts ; ks++ ){
00319          if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
00320             initialize = 1 ;
00321             fprintf(stderr,"Inadequate time series #%d in L1F plugin\n\a",ks+1) ;
00322             return ;
00323          }
00324          tsar = MRI_FLOAT_PTR(tsim[ks]) ;
00325          for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = tsar[ii+ignore] ;
00326          ir++ ;
00327       }
00328 
00329       /* make space for fit vector */
00330 
00331       fit = (float *) malloc(sizeof(float)*nref) ;
00332 
00333       initialize = 0 ;
00334    }
00335 
00336    /** find L1 fit coefficients **/
00337 
00338    cls = cl1_solve( nlen , nref , vec+ignore , ref , fit,0 ) ;
00339 
00340    if( cls < 0.0 ) return ;  /* bad fit */
00341 
00342    for( ii=0 ; ii < nlen ; ii++ ){
00343       val = 0.0 ;
00344       for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
00345 
00346       vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
00347    }
00348 
00349    /** 22 Apr 1997: create label if desired by AFNI         **/
00350    /** [This is in static storage, since AFNI will copy it] **/
00351 
00352    if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00353 
00354       lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00355 
00356       /** for each reference, make a string into sbuf **/
00357 
00358       ir = 0 ;
00359       sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
00360       strcat(lbuf,sbuf) ;
00361 
00362       for( ; ir <= polort ; ){
00363          sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir++] ) ;
00364          strcat(lbuf,sbuf) ;
00365       }
00366 
00367       for( ks=0 ; ks < nrsin ; ks++ ){
00368          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00369             fac = sinper[ks] / (jh+1) ;
00370             sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00371             strcat(lbuf,sbuf) ;
00372             sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00373             strcat(lbuf,sbuf) ;
00374          }
00375       }
00376 
00377       for( ks=0 ; ks < nrts ; ks++ ){
00378          sprintf(sbuf,"Coef of %s = %g\n" , tsim[ks]->name , fit[ir++] ) ;
00379          strcat(lbuf,sbuf) ;
00380       }
00381 
00382       *label = lbuf ;  /* send address of lbuf back in what label points to */
00383    }
00384 
00385    return ;
00386 }