00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #define PROGRAM_NAME "plug_wavelets"
00026 #define PROGRAM_AUTHOR "B. Douglas Ward"
00027 #define PROGRAM_INITIAL "28 March 2000"
00028 #define PROGRAM_LATEST "02 December 2002"
00029
00030
00031
00032
00033
00034
00035 #include "afni.h"
00036 #include "Wavelets.h"
00037 #include "Wavelets.c"
00038
00039
00040
00041
00042
00043
00044 #define MAX_NAME_LENGTH THD_MAX_NAME
00045 #define MAX_FILTERS 20
00046 #define MAX_BAND 20
00047
00048
00049
00050 #ifndef ALLOW_PLUGINS
00051 # error "Plugins not properly set up -- see machdep.h"
00052 #endif
00053
00054
00055
00056
00057 static char helpstring[] =
00058 " Purpose: Wavelet Analysis of FMRI time series data. \n"
00059 " \n"
00060 " Control: Wavelet = Type of wavelet to be used in analysis \n"
00061 " NFirst = Number of first time series point to use. \n"
00062 " NLast = Number of last time series point to use. \n"
00063 " \n"
00064 " Filter: Type = Type of filtering to apply in this \n"
00065 " time-frequency window: \n"
00066 " Stop = set wavelet coefficients to zero \n"
00067 " Baseline = assign wavelet coefficients to baseline model \n"
00068 " Signal = assign wavelet coefficients to signal model \n"
00069 " \n"
00070 " Band = Frequency band at which to apply filtering. \n"
00071 " Min TR = Minimum value for time window (in TR) \n"
00072 " Max TR = Maximum value for time window (in TR) \n"
00073 " \n"
00074 "For more information, see 'Wavelet Analysis of FMRI Time Series Data' \n"
00075 "which is contained in file 3dWavelets.ps of the AFNI distribution. \n"
00076 ;
00077
00078
00079
00080 #define NFILTER 3
00081 static char * filter_strings[NFILTER] = {"Stop", "Baseline", "Signal"};
00082
00083
00084
00085
00086 char * WA_main( PLUGIN_interface * ) ;
00087
00088 void WA_fwt (int nt, double to, double dt, float * vec, char ** label);
00089 void WA_fit (int nt, double to, double dt, float * vec, char ** label);
00090 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label);
00091 void WA_err (int nt, double to, double dt, float * vec, char ** label);
00092
00093
00094
00095
00096
00097 static PLUGIN_interface * global_plint = NULL;
00098
00099 static int plug_wavelet_type=0;
00100 static int plug_NFirst=0;
00101 static int plug_NLast=32767;
00102 static int plug_initialize=0;
00103 static int plug_prev_nt=0;
00104 static int plug_filter_type=0;
00105
00106 static int plug_num_stop_filters = 0;
00107 static int plug_stop_band[MAX_FILTERS];
00108 static int plug_stop_mintr[MAX_FILTERS];
00109 static int plug_stop_maxtr[MAX_FILTERS];
00110 static float * plug_stop_filter = NULL;
00111
00112 static int plug_num_base_filters = 0;
00113 static int plug_base_band[MAX_FILTERS];
00114 static int plug_base_mintr[MAX_FILTERS];
00115 static int plug_base_maxtr[MAX_FILTERS];
00116 static float * plug_base_filter = NULL;
00117
00118
00119 static int plug_num_sgnl_filters = 0;
00120 static int plug_sgnl_band[MAX_FILTERS];
00121 static int plug_sgnl_mintr[MAX_FILTERS];
00122 static int plug_sgnl_maxtr[MAX_FILTERS];
00123 static float * plug_sgnl_filter = NULL;
00124
00125
00126
00127
00128
00129
00130
00131
00132 DEFINE_PLUGIN_PROTOTYPE
00133
00134 PLUGIN_interface * PLUGIN_init( int ncall )
00135 {
00136 PLUGIN_interface * plint ;
00137 int is;
00138
00139
00140 if( ncall > 0 ) return NULL ;
00141
00142
00143
00144
00145
00146
00147 plint = PLUTO_new_interface ("Wavelets" ,
00148 "Wavelet Analysis of Time Series Data" ,
00149 helpstring, PLUGIN_CALL_VIA_MENU, WA_main);
00150
00151 global_plint = plint ;
00152
00153 PLUTO_add_hint (plint,
00154 "Control Wavelet Analysis Functions");
00155
00156 PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00157
00158
00159
00160 for (is =0; is < MAX_FILTERS; is++)
00161 {
00162 plug_stop_band[is] = 0;
00163 plug_stop_mintr[is] = 0.0;
00164 plug_stop_maxtr[is] = 0.0;
00165 plug_base_band[is] = 0;
00166 plug_base_mintr[is] = 0.0;
00167 plug_base_maxtr[is] = 0.0;
00168 plug_sgnl_band[is] = 0;
00169 plug_sgnl_mintr[is] = 0.0;
00170 plug_sgnl_maxtr[is] = 0.0;
00171 }
00172
00173
00174
00175 PLUTO_add_option (plint, "Control", "Control", TRUE);
00176 PLUTO_add_string (plint, "Wavelet", MAX_WAVELET_TYPE, WAVELET_TYPE_name,
00177 plug_wavelet_type);
00178 PLUTO_add_number (plint, "NFirst", 0, 32767, 0, plug_NFirst, TRUE);
00179 PLUTO_add_number (plint, "NLast", 0, 32767, 0, plug_NLast, TRUE);
00180
00181
00182
00183 for (is = 0; is < MAX_FILTERS; is++)
00184 {
00185 PLUTO_add_option (plint, "Filter", "Filter", FALSE);
00186 PLUTO_add_string (plint, "Type", NFILTER, filter_strings,
00187 plug_filter_type);
00188 PLUTO_add_number (plint, "Band", -1, MAX_BAND, 0, 0, TRUE);
00189 PLUTO_add_number (plint, "Min TR", 0, 10000, 0, 0, TRUE);
00190 PLUTO_add_number (plint, "Max TR", 0, 10000, 0, 0, TRUE);
00191 }
00192
00193
00194
00195 PLUTO_register_1D_funcstr ("WA_FWT", WA_fwt);
00196 PLUTO_register_1D_funcstr ("WA_Fit", WA_fit);
00197 PLUTO_register_1D_funcstr ("WA_Sgnl", WA_sgnl);
00198 PLUTO_register_1D_funcstr ("WA_Err", WA_err);
00199
00200
00201 return plint ;
00202 }
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 char * WA_main( PLUGIN_interface * plint )
00213 {
00214 char * str;
00215 int is;
00216
00217
00218
00219 plug_initialize = 0;
00220
00221
00222
00223 PLUTO_next_option (plint);
00224 str = PLUTO_get_string (plint);
00225 plug_wavelet_type = PLUTO_string_index (str, MAX_WAVELET_TYPE,
00226 WAVELET_TYPE_name);
00227 plug_NFirst = PLUTO_get_number (plint);
00228 plug_NLast = PLUTO_get_number (plint);
00229
00230
00231
00232 plug_num_stop_filters = 0;
00233 plug_num_base_filters = 0;
00234 plug_num_sgnl_filters = 0;
00235
00236 do
00237 {
00238 str = PLUTO_get_optiontag(plint);
00239 if (str == NULL) break;
00240 if (strcmp (str, "Filter") != 0)
00241 return "************************\n"
00242 "Illegal optiontag found!\n"
00243 "************************";
00244
00245
00246
00247 str = PLUTO_get_string (plint);
00248 plug_filter_type = PLUTO_string_index (str, NFILTER, filter_strings);
00249
00250 switch (plug_filter_type)
00251 {
00252 case 0:
00253 {
00254 plug_stop_band[plug_num_stop_filters] = PLUTO_get_number(plint);
00255 plug_stop_mintr[plug_num_stop_filters]
00256 = PLUTO_get_number(plint);
00257 plug_stop_maxtr[plug_num_stop_filters]
00258 = PLUTO_get_number(plint);
00259
00260 if (plug_stop_mintr[plug_num_stop_filters] >
00261 plug_stop_maxtr[plug_num_stop_filters])
00262 return "*************************\n"
00263 "Require Min TR <= Max TR \n"
00264 "*************************" ;
00265
00266 plug_num_stop_filters++;
00267 break;
00268 }
00269
00270 case 1:
00271 {
00272 plug_base_band[plug_num_base_filters] = PLUTO_get_number(plint);
00273 plug_base_mintr[plug_num_base_filters]
00274 = PLUTO_get_number(plint);
00275 plug_base_maxtr[plug_num_base_filters]
00276 = PLUTO_get_number(plint);
00277
00278 if (plug_base_mintr[plug_num_base_filters] >
00279 plug_base_maxtr[plug_num_base_filters])
00280 return "*************************\n"
00281 "Require Min TR <= Max TR \n"
00282 "*************************" ;
00283
00284 plug_num_base_filters++;
00285 break;
00286 }
00287
00288 case 2:
00289 {
00290 plug_sgnl_band[plug_num_sgnl_filters]=PLUTO_get_number(plint);
00291 plug_sgnl_mintr[plug_num_sgnl_filters]
00292 = PLUTO_get_number(plint);
00293 plug_sgnl_maxtr[plug_num_sgnl_filters]
00294 = PLUTO_get_number(plint);
00295
00296 if (plug_sgnl_mintr[plug_num_sgnl_filters] >
00297 plug_sgnl_maxtr[plug_num_sgnl_filters])
00298 return "*************************\n"
00299 "Require Min TR <= Max TR \n"
00300 "*************************" ;
00301
00302 plug_num_sgnl_filters++;
00303 break;
00304 }
00305
00306 }
00307 }
00308 while (1);
00309
00310
00311
00312 printf ("\n\n");
00313 printf ("Program: %s \n", PROGRAM_NAME);
00314 printf ("Author: %s \n", PROGRAM_AUTHOR);
00315 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00316 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00317 printf ("\n");
00318
00319
00320 printf ("\nControls: \n");
00321 printf ("Wavelet Type = %10s \n", WAVELET_TYPE_name[plug_wavelet_type]);
00322 printf ("NFirst = %10d \n", plug_NFirst);
00323 printf ("NLast = %10d \n", plug_NLast);
00324
00325 for (is = 0; is < plug_num_stop_filters; is++)
00326 {
00327 printf ("\nStop Filter: Band = %4d ", plug_stop_band[is]);
00328 printf ("Min. TR = %4d Max. TR = %4d \n",
00329 plug_stop_mintr[is], plug_stop_maxtr[is]);
00330 }
00331
00332 for (is = 0; is < plug_num_base_filters; is++)
00333 {
00334 printf ("\nBaseline Filter: Band = %4d ", plug_base_band[is]);
00335 printf ("Min. TR = %4d Max. TR = %4d \n",
00336 plug_base_mintr[is], plug_base_maxtr[is]);
00337 }
00338
00339 for (is = 0; is < plug_num_sgnl_filters; is++)
00340 {
00341 printf ("\nSignal Filter: Band = %4d ", plug_sgnl_band[is]);
00342 printf ("Min. TR = %4d Max. TR = %4d \n",
00343 plug_sgnl_mintr[is], plug_sgnl_maxtr[is]);
00344 }
00345
00346
00347
00348 plug_initialize = 1 ;
00349 plug_prev_nt = 0;
00350
00351 return NULL ;
00352 }
00353
00354
00355
00356
00357
00358
00359
00360 int calculate_results
00361 (
00362 int nt,
00363 float * vec,
00364 int * NFirst,
00365 int * NLast,
00366 char ** label,
00367 float ** coefts,
00368 float ** fitts,
00369 float ** sgnlts,
00370 float ** errts
00371 )
00372
00373 {
00374 float * ts_array;
00375 float * coef;
00376 float sse_base;
00377 float sse_full;
00378 float ffull;
00379 float rfull;
00380
00381 int N;
00382 int f;
00383 int p;
00384 int q;
00385 int n;
00386 int i;
00387 int ok = 1;
00388
00389
00390
00391 if (! plug_initialize) return (0);
00392
00393
00394
00395 *NFirst = plug_NFirst;
00396
00397 *NLast = plug_NLast;
00398 if (*NLast > nt-1) *NLast = nt-1;
00399
00400 N = *NLast - *NFirst + 1;
00401 N = powerof2(my_log2(N));
00402 *NLast = N + *NFirst - 1;
00403
00404
00405 plug_stop_filter =
00406 FWT_1d_stop_filter (plug_num_stop_filters, plug_stop_band,
00407 plug_stop_mintr, plug_stop_maxtr, *NFirst, N);
00408
00409 plug_base_filter =
00410 FWT_1d_pass_filter (plug_num_base_filters, plug_base_band,
00411 plug_base_mintr, plug_base_maxtr, *NFirst, N);
00412
00413 plug_sgnl_filter =
00414 FWT_1d_pass_filter (plug_num_sgnl_filters, plug_sgnl_band,
00415 plug_sgnl_mintr, plug_sgnl_maxtr, *NFirst, N);
00416
00417 f = 0;
00418 for (i = 0; i < N; i++)
00419 if (plug_stop_filter[i] == 0.0)
00420 {
00421 f++;
00422 plug_base_filter[i] = 0.0;
00423 plug_sgnl_filter[i] = 0.0;
00424 }
00425
00426 q = 0;
00427 for (i = 0; i < N; i++)
00428 if (plug_base_filter[i] == 1.0)
00429 {
00430 q++;
00431 plug_sgnl_filter[i] = 1.0;
00432 }
00433
00434 p = 0;
00435 for (i = 0; i < N; i++)
00436 if (plug_sgnl_filter[i] == 1.0)
00437 {
00438 p++;
00439 }
00440
00441
00442
00443 coef = (float *) malloc (sizeof(float) * p); MTEST (coef);
00444 *coefts = (float *) malloc (sizeof(float) * N); MTEST (coefts);
00445 *fitts = (float *) malloc (sizeof(float) * N); MTEST (fitts);
00446 *sgnlts = (float *) malloc (sizeof(float) * N); MTEST (sgnlts);
00447 *errts = (float *) malloc (sizeof(float) * N); MTEST (errts);
00448
00449
00450
00451 ts_array = vec + *NFirst;
00452
00453
00454
00455 wavelet_analysis (plug_wavelet_type,
00456 f, plug_stop_filter,
00457 q, plug_base_filter,
00458 p, plug_sgnl_filter,
00459 N, ts_array, coef,
00460 &sse_base, &sse_full, &ffull, &rfull,
00461 *coefts, *fitts, *sgnlts, *errts);
00462
00463
00464
00465 printf ("\nResults for Voxel: \n");
00466 report_results (N, *NFirst, f, p, q,
00467 plug_base_filter, plug_sgnl_filter,
00468 coef, sse_base, sse_full, ffull, rfull,
00469 label);
00470 printf ("%s \n", *label);
00471
00472 plug_prev_nt = nt;
00473
00474
00475
00476 free (plug_stop_filter); plug_stop_filter = NULL;
00477 free (plug_base_filter); plug_base_filter = NULL;
00478 free (plug_sgnl_filter); plug_sgnl_filter = NULL;
00479 free (coef); coef = NULL;
00480
00481
00482 return (ok);
00483 }
00484
00485
00486
00487
00488
00489
00490
00491 void WA_fwt (int nt, double to, double dt, float * vec, char ** label)
00492
00493 {
00494 int NFirst;
00495 int NLast;
00496 int n;
00497 int ok;
00498 float * coefts = NULL;
00499 float * fitts = NULL;
00500 float * sgnlts = NULL;
00501 float * errts = NULL;
00502
00503
00504
00505 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
00506 &coefts, &fitts, &sgnlts, &errts);
00507 if (!ok)
00508 {
00509 for (n = 0; n < nt; n++) vec[n] = 0.0;
00510 return;
00511 }
00512
00513
00514
00515 for (n = NFirst; n <= NLast; n++)
00516 {
00517 vec[n] = coefts[n-NFirst];
00518 }
00519
00520 for (n = 0; n < NFirst; n++)
00521 vec[n] = 0.0;
00522
00523 for (n = NLast+1; n < nt; n++)
00524 vec[n] = 0.0;
00525
00526
00527
00528 free (coefts); coefts = NULL;
00529 free (fitts); fitts = NULL;
00530 free (sgnlts); sgnlts = NULL;
00531 free (errts); errts = NULL;
00532
00533
00534 return;
00535 }
00536
00537
00538
00539
00540
00541
00542
00543 void WA_fit (int nt, double to, double dt, float * vec, char ** label)
00544
00545 {
00546 int NFirst;
00547 int NLast;
00548 int n;
00549 int ok;
00550 float * coefts = NULL;
00551 float * fitts = NULL;
00552 float * sgnlts = NULL;
00553 float * errts = NULL;
00554
00555
00556
00557 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
00558 &coefts, &fitts, &sgnlts, &errts);
00559 if (!ok)
00560 {
00561 for (n = 0; n < nt; n++) vec[n] = 0.0;
00562 return;
00563 }
00564
00565
00566
00567 for (n = NFirst; n <= NLast; n++)
00568 {
00569 vec[n] = fitts[n-NFirst];
00570 }
00571
00572 for (n = 0; n < NFirst; n++)
00573 vec[n] = vec[NFirst];
00574
00575 for (n = NLast+1; n < nt; n++)
00576 vec[n] = vec[NLast];
00577
00578
00579
00580 free (coefts); coefts = NULL;
00581 free (fitts); fitts = NULL;
00582 free (sgnlts); sgnlts = NULL;
00583 free (errts); errts = NULL;
00584
00585
00586 return;
00587 }
00588
00589
00590
00591
00592
00593
00594
00595 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label)
00596
00597 {
00598 int NFirst;
00599 int NLast;
00600 int n;
00601 int ok;
00602 float * coefts = NULL;
00603 float * fitts = NULL;
00604 float * sgnlts = NULL;
00605 float * errts = NULL;
00606
00607
00608
00609 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
00610 &coefts, &fitts, &sgnlts, &errts);
00611 if (!ok)
00612 {
00613 for (n = 0; n < nt; n++) vec[n] = 0.0;
00614 return;
00615 }
00616
00617
00618
00619 for (n = NFirst; n <= NLast; n++)
00620 {
00621 vec[n] = sgnlts[n-NFirst];
00622 }
00623
00624 for (n = 0; n < NFirst; n++)
00625 vec[n] = 0.0;
00626
00627 for (n = NLast+1; n < nt; n++)
00628 vec[n] = 0.0;
00629
00630
00631
00632 free (coefts); coefts = NULL;
00633 free (fitts); fitts = NULL;
00634 free (sgnlts); sgnlts = NULL;
00635 free (errts); errts = NULL;
00636
00637
00638 return;
00639 }
00640
00641
00642
00643
00644
00645
00646
00647 void WA_err (int nt, double to, double dt, float * vec, char ** label)
00648
00649 {
00650 int NFirst;
00651 int NLast;
00652 int n;
00653 int ok;
00654 float * coefts = NULL;
00655 float * fitts = NULL;
00656 float * sgnlts = NULL;
00657 float * errts = NULL;
00658
00659
00660
00661 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
00662 &coefts, &fitts, &sgnlts, &errts);
00663 if (!ok)
00664 {
00665 for (n = 0; n < nt; n++) vec[n] = 0.0;
00666 return;
00667 }
00668
00669
00670
00671 for (n = NFirst; n <= NLast; n++)
00672 {
00673 vec[n] = errts[n-NFirst];
00674 }
00675
00676 for (n = 0; n < NFirst; n++)
00677 vec[n] = 0.0;
00678
00679 for (n = NLast+1; n < nt; n++)
00680 vec[n] = 0.0;
00681
00682
00683
00684 free (coefts); coefts = NULL;
00685 free (fitts); fitts = NULL;
00686 free (sgnlts); sgnlts = NULL;
00687 free (errts); errts = NULL;
00688
00689
00690 return;
00691 }
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708