00001 #include "afni.h"
00002
00003 #ifndef ALLOW_PLUGINS
00004 # error "Plugins not properly set up -- see machdep.h"
00005 #endif
00006
00007
00008
00009
00010
00011
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
00028
00029 #define NBASE 4
00030 static char * baseline_strings[NBASE] = { "Constant", "Linear", "Quadratic", "Cubic" } ;
00031
00032
00033
00034 char * L1F_main( PLUGIN_interface * ) ;
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
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
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 DEFINE_PLUGIN_PROTOTYPE
00069
00070 PLUGIN_interface * PLUGIN_init( int ncall )
00071 {
00072 int ii ;
00073 PLUGIN_interface * plint ;
00074
00075 if( ncall > 0 ) return NULL ;
00076
00077
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 ;
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" ) ;
00091
00092
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
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
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
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
00125
00126
00127
00128
00129 char * L1F_main( PLUGIN_interface * plint )
00130 {
00131 char * str ;
00132 int ii ;
00133 float * tsar ;
00134
00135
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
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++ ) ;
00172 ignore = ii ;
00173 nrts++ ;
00174
00175 } else {
00176 return "************************\n"
00177 "Illegal optiontag found!\n"
00178 "************************" ;
00179 }
00180 } while(1) ;
00181
00182
00183
00184 initialize = 1 ;
00185
00186
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
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] ;
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
00233
00234 nref = (polort+1) + nrts ;
00235 for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00236
00237
00238
00239 nlen = nt - ignore ;
00240
00241 if( nlen <= nref ) return ;
00242
00243
00244
00245
00246
00247
00248
00249 if( nlen != nlen_old || nref != nref_old ||
00250 initialize || (dt != dt_old && nrsin > 0) ){
00251
00252
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
00261
00262 ref = (float **) malloc( sizeof(float *) * nref ) ;
00263 if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00264 return;
00265
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
00273 }
00274 }
00275 nlen_old = nlen ;
00276 nref_old = nref ;
00277 dt_old = dt ;
00278
00279
00280
00281
00282
00283 for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
00284
00285 ir = 1 ;
00286 if( polort > 0 ){
00287
00288
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
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
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
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
00330
00331 fit = (float *) malloc(sizeof(float)*nref) ;
00332
00333 initialize = 0 ;
00334 }
00335
00336
00337
00338 cls = cl1_solve( nlen , nref , vec+ignore , ref , fit,0 ) ;
00339
00340 if( cls < 0.0 ) return ;
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
00350
00351
00352 if( label != NULL ){
00353
00354 lbuf[0] = '\0' ;
00355
00356
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 ;
00383 }
00384
00385 return ;
00386 }