1 #include "mrilib.h"
2 #include "matrix.h"
3 #include "plug_delay_V2.h"
4
5 static void RP_tsfunc( double tzero, double tdelta ,
6 int npts, float ts[],
7 double ts_mean , double ts_slope ,
8 void *ud, int nbriks, float *val );
9
10 static void DEL_tsfunc( double tzero, double tdelta ,
11 int npts, float ts[],
12 double ts_mean , double ts_slope ,
13 void *ud, int nbriks, float *val );
14
15 typedef enum { FFT_PHASE = 1, HILB_DELAY = 2 } PHASE_METHODS;
16
17 typedef enum {
18 NEG=-3, CONT=-2, CCW=-1,
19 NOT_SET=0,
20 CLW=1, EXP=2, POS=3
21 } PHAZE_DIRS;
22 typedef enum {
23 k_EXP=0, k_CON, k_CLW, k_CCW, k_N } STIM_TYPES;
24 typedef enum {
25 ECC=0, POL=1 } FIELD_PARAMS;
26
27 typedef struct {
28 int pmeth; /* 1 --> phase, 2 -->delay */
29 /* to set outside of processor function */
30 float ftap; /* taper fraction at ends of data */
31 float dt; /* TR */
32 float pre; /* period in seconds, before stimulus began */
33 float Fstim[2]; /* stimulus fundamental freq. in Hz */
34 int nvals; /* number of points in time series */
35 int verb;
36 THD_3dim_dataset *iset;
37 int *vox; /* vox[k] is the 1D index of the kth voxel in mask*/
38 int nmask; /* total number of voxels in mask */
39 char *prefix;
40 char *oext; /* output extension */
41 int dir; /* +1 for CLW, -1 for CCW
42 +2 for Exp. -2 for contracting*/
43 int n[2]; /* number of rings, wedges */
44 int spectra;
45 int fixsum;
46 int ort_adj; /* number of orts removed from data */
47
48 /* to be set upon first calling of processor function */
49 int nfft; /* number of points in fft */
50 float fstep; /* Delta freq. Hz*/
51 float Fresp[2]; /* FMRI response fundamental freq. in Hz */
52 int stk[2]; /* floor, and ceil, of index of fresp in fft array */
53 float stw[2]; /* weight of floor and ceil for index of fresp in fft array */
54 THD_3dim_dataset *phz;
55 THD_3dim_dataset *amp;
56 int dof[2];
57
58 /* ops specific to hilbert phase estimation */
59 float **rvec; /* reference time series */
60 int rvec_len; /* number of time points in rvec */
61 int rvec_num; /* number of rvecs */
62 int iref; /* which rvec are we dealing with? */
63 int Dsamp; /* correct slice timing offset */
64
65 } RP_UD;
66
67 #define PHASE_R2D(a) ( (a)*180.0/PI )
68 #define PHASE_360(a) ( (a) < 0.0 ? 360.0+(a):(a) )
Phase_Type_to_Dir(int p)69 int Phase_Type_to_Dir(int p) {
70 switch(p) {
71 case k_CON:
72 return(CONT);
73 case k_EXP:
74 return(EXP);
75 case k_CLW:
76 return(CLW);
77 case k_CCW:
78 return(CCW);
79 default:
80 return(NOT_SET);
81 }
82 }
83
Phase_Dir_to_Type(int p)84 int Phase_Dir_to_Type(int p) {
85 switch(p) {
86 case CONT:
87 return(k_CON);
88 case EXP:
89 return(k_EXP);
90 case CLW:
91 return(k_CLW);
92 case CCW:
93 return(k_CCW);
94 default:
95 return(-1);
96 }
97 }
98
Phase_Method_string(int p)99 char * Phase_Method_string(int p) {
100 static char ps[64]={"ERROR"};
101 switch(p) {
102 case FFT_PHASE:
103 sprintf(ps,"FFT-Phase");
104 break;
105 case HILB_DELAY:
106 sprintf(ps,"Hilbert-delay");
107 break;
108 default:
109 sprintf(ps,"Calamity");
110 break;
111 }
112 return(ps);
113 }
114
Phase_Dirs_string(int p)115 char * Phase_Dirs_string(int p) {
116 static char ps[64]={"ERROR"};
117 switch(p) {
118 case CONT:
119 sprintf(ps,"Contracting (Eccentricity -)");
120 break;
121 case CCW:
122 sprintf(ps,"Counter clockwise (Polar -)");
123 break;
124 case CLW:
125 sprintf(ps,"Clockwise (Polar +)");
126 break;
127 case EXP:
128 sprintf(ps,"Expanding (Eccentricity +)");
129 break;
130 case POS:
131 sprintf(ps,"Positive + (Polar, or Ecc)");
132 break;
133 case NEG:
134 sprintf(ps,"Negative - (Polar, or Ecc)");
135 break;
136 case NOT_SET:
137 sprintf(ps,"Not_set");
138 break;
139 default:
140 sprintf(ps,"Horreur");
141 break;
142 }
143 return(ps);
144 }
Phase_Dirs_lbl(int p)145 char * Phase_Dirs_lbl(int p) {
146 static char ps[64]={"ERROR"};
147 switch(p) {
148 case CONT:
149 sprintf(ps,"ecc-");
150 break;
151 case CCW:
152 sprintf(ps,"pol-");
153 break;
154 case CLW:
155 sprintf(ps,"pol+");
156 break;
157 case EXP:
158 sprintf(ps,"ecc+");
159 break;
160 case POS:
161 sprintf(ps,"+");
162 break;
163 case NEG:
164 sprintf(ps,"-");
165 break;
166 case NOT_SET:
167 sprintf(ps,".");
168 break;
169 default:
170 sprintf(ps,"Horreur");
171 break;
172 }
173 return(ps);
174 }
Phase_Dirs_ulbl(int p)175 char * Phase_Dirs_ulbl(int p) {
176 static char ps[64]={"ERROR"};
177 switch(p) {
178 case CONT:
179 case EXP:
180 sprintf(ps,"ecc");
181 break;
182 case CCW:
183 case CLW:
184 sprintf(ps,"pol");
185 break;
186 default:
187 sprintf(ps,"Horreur");
188 break;
189 }
190 return(ps);
191 }
192
Dir_is_eccentricity(d)193 int Dir_is_eccentricity(d) {
194 if (d == CONT || d == EXP) return(1);
195 return(0);
196 }
197
Dir_is_polar(d)198 int Dir_is_polar(d) {
199 if (d == CLW || d == CCW) return(1);
200 return(0);
201 }
202
203
Dir2Type(p)204 int Dir2Type(p) {
205 switch (p){
206 case CONT:
207 case EXP:
208 return(ECC);
209 break;
210 case CCW:
211 case CLW:
212 return(POL);
213 break;
214 default:
215 return(-1);
216 break;
217 }
218 }
219
Show_RP_UD(RP_UD * u,char * str)220 Show_RP_UD(RP_UD *u, char *str) {
221 if (str) {
222 fprintf(stderr,"%s", str);
223 }
224 fprintf(stderr,"Phase Method %d, %s\n"
225 " nfft=%d, fstep=%.3f\n"
226 " nvals=%d\n"
227 " stk=[%d,%d]\n"
228 " stw=[%.3f, %.3f]\n"
229 " ftap=%f\n"
230 " dt=%f, pre = %f\n"
231 " Eccentricity: fstim=%f, fresp=%f\n"
232 " Polar: fstim=%f, fresp=%f\n"
233 " dof=[%d %d], ort_adj=%d\n"
234 " verb=%d, spectra=%d\n"
235 " iset=%p\n"
236 " vox=%p, nmask=%d\n"
237 " phz=%p\n"
238 " amp=%p\n"
239 " prefix=%s, oext=%s\n"
240 " %d wedges, %d rings\n"
241 " fixsum=%d\n"
242 " This call dir=%d (%s) \n"
243 " Dsamp = %d, rvec=%p, %d vals, %d refs, iref = %d \n"
244 ,
245 u->pmeth, Phase_Method_string(u->pmeth),
246 u->nfft, u->fstep, u->nvals,
247 u->stk[0], u->stk[1],
248 u->stw[0], u->stw[1],
249 u->ftap, u->dt, u->pre,
250 u->Fstim[ECC], u->Fresp[ECC],
251 u->Fstim[POL], u->Fresp[POL],
252 u->dof[0], u->dof[1], u->ort_adj,
253 u->verb, u->spectra,
254 u->iset, u->vox, u->nmask,
255 u->amp, u->phz,
256 u->prefix ? u->prefix:"NULL",
257 u->oext ? u->oext:"NULL",
258 u->n[POL], u->n[ECC], u->fixsum,
259 u->dir, Phase_Dirs_string(u->dir),
260 u->Dsamp, u->rvec, u->rvec_len, u->rvec_num, u->iref);
261 }
262
SetFreqBin(float fresp,float fstep,int stk[],float stw[],int nfft)263 int SetFreqBin(float fresp, float fstep, int stk[], float stw[], int nfft) {
264 float stf = fresp/fstep; /* stimulus freq. index */
265
266 stk[0] = (int)floor(stf); /* floor of freq. index */
267 stw[0] = (1.0 - stf + stk[0]); /*floor weight*/
268 stk[1] = (int)ceil (stf); /* ceil of freq. index */
269 stw[1] = (1.0 - stk[1] + stf ); /*ceil weight*/
270 if (stk[0] > nfft/2 || stk[1] > nfft/2) {
271 ERROR_message("Frequency indices %d and/or %d outside max of %d\n",
272 stk[0], stk[1], nfft/2);
273 return(0);
274 }
275 return(1);
276 }
277
278
Combine_Opposites(THD_3dim_dataset * dset1,THD_3dim_dataset * dset2,RP_UD * rpud)279 THD_3dim_dataset * Combine_Opposites(THD_3dim_dataset *dset1,
280 THD_3dim_dataset *dset2,
281 RP_UD *rpud)
282 {
283 float *phi1=NULL, *phi2=NULL, *xc1=NULL, *xc2=NULL, *mxxc=NULL,
284 *sum=NULL, *dif=NULL;
285 int nvox, i;
286 THD_3dim_dataset *oset = NULL;
287 char stmp[256+strlen(rpud->prefix)];
288 double radpersec = 0.0, n=1.0;
289
290 radpersec = (360.0*rpud->Fresp[Dir2Type(rpud->dir)]);
291
292 n = (float)rpud->n[Dir2Type(rpud->dir)];
293
294 nvox = DSET_NVOX(dset1);
295 phi1 = (float *)calloc(nvox, sizeof(float));
296 xc1 = (float *)calloc(nvox, sizeof(float));
297 phi2 = (float *)calloc(nvox, sizeof(float));
298 xc2 = (float *)calloc(nvox, sizeof(float));
299 sum = (float *)calloc(nvox, sizeof(float));
300 dif = (float *)calloc(nvox, sizeof(float));
301 mxxc = (float *)calloc(nvox, sizeof(float));
302
303 EDIT_coerce_scale_type( nvox , DSET_BRICK_FACTOR(dset1,0) ,
304 DSET_BRICK_TYPE(dset1,0),
305 DSET_ARRAY(dset1, 0) , /* input */
306 MRI_float, phi1 ) ;
307 EDIT_coerce_scale_type( nvox , DSET_BRICK_FACTOR(dset1,1) ,
308 DSET_BRICK_TYPE(dset1,1),
309 DSET_ARRAY(dset1, 1) , /* input */
310 MRI_float, xc1 ) ;
311 EDIT_coerce_scale_type( nvox , DSET_BRICK_FACTOR(dset2,0) ,
312 DSET_BRICK_TYPE(dset2,0),
313 DSET_ARRAY(dset2, 0) , /* input */
314 MRI_float, phi2 ) ;
315 EDIT_coerce_scale_type( nvox , DSET_BRICK_FACTOR(dset2,1) ,
316 DSET_BRICK_TYPE(dset2,1),
317 DSET_ARRAY(dset2, 1) , /* input */
318 MRI_float, xc2 ) ;
319
320 for (i=0; i<nvox;++i) {
321 dif[i] = (phi1[i]-phi2[i])/2.0;
322 if (ABS(dif[i]) < 90/n || !rpud->fixsum) { /* should this be 90 / n ? */
323 sum[i] = (phi1[i]+phi2[i])/2.0;
324 } else { /* Too big a difference in phase, likely
325 summing about zero degrees where wrapping from
326 360 occurs*/
327 sum[i] = (phi1[i]+phi2[i])/2.0-180.0/n;
328 if (sum[i]<0) sum[i] = 360/n + sum[i];
329 }
330 /* now change difference of phase from degrees to seconds for output */
331 /* dif[i] /= radpersec; */
332 /* keep track of max coef == union masking */
333 if (xc2[i] > xc1[i]) mxxc[i] = xc2[i]; else mxxc[i] = xc1[i];
334 }
335
336 /* put the results in a dset for output */
337 oset = EDIT_empty_copy( dset1 ) ;
338 sprintf(stmp,"%s.%s.field%s",
339 rpud->prefix, Phase_Dirs_ulbl(rpud->dir),rpud->oext);
340 EDIT_dset_items( oset ,
341 ADN_prefix , stmp,
342 ADN_datum_all, MRI_float ,
343 ADN_nvals , 3 ,
344 ADN_none ) ;
345
346 EDIT_substitute_brick( oset , 0 , MRI_float , sum ) ; /* do not free sum */
347 EDIT_substitute_brick( oset , 1 , MRI_float , dif ) ; /* do not free dif */
348 EDIT_substitute_brick( oset , 2 , MRI_float , mxxc ) ;/* do not free mxxc */
349 if (Dir_is_eccentricity(rpud->dir)) {
350 EDIT_BRICK_LABEL(oset , 0, "Eccentricity");
351 } else if (Dir_is_polar(rpud->dir)) {
352 EDIT_BRICK_LABEL(oset , 0, "Polar Angle");
353 } else {
354 ERROR_message("rpud->dir makes no sense here");
355 }
356 EDIT_BRICK_LABEL(oset , 1, "Hemo. Offset");
357 if (rpud->pmeth == FFT_PHASE) {
358 sprintf(stmp,"Max.PwR@%.3fHz", rpud->Fresp[Dir2Type(rpud->dir)]);
359 EDIT_BRICK_LABEL(oset , 2, stmp);
360 } else {
361 EDIT_BRICK_LABEL(oset , 2, "Max.Corr.Coef.");
362 /* Could assume dset1 and dset2 have the same length series, and dofs
363 and do:
364 EDIT_BRICK_TO_FICO(dset1, 2, rpud->rvec_len,
365 rpud->dof[0],rpud->dof[1]);
366 But it is safer to deprive users of significance levels which
367 don't quite apply with the max operation */
368 }
369 free(phi1); phi1 = NULL;
370 free(phi2); phi2 = NULL;
371 free(xc1); xc1 = NULL;
372 free(xc2); xc2 = NULL;
373
374 return(oset);
375 }
376
RP_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)377 static void RP_tsfunc( double tzero, double tdelta ,
378 int npts, float ts[],
379 double ts_mean , double ts_slope ,
380 void *ud, int nbriks, float *val )
381 {
382 static int nvox , ncall , N_nzfreq=0, *stimharm=NULL;
383 static byte *nzfreq=NULL;
384 static float *xtap=NULL, *mag=NULL, *phz=NULL, stf;
385 static complex *comp_array=NULL;
386 int ii , jj;
387 double nzpow=0.0, preoff = 0.0;
388 RP_UD *rpud = (RP_UD *)ud;
389 static int st = -1;
390 /** is this a "notification"? **/
391
392 if (!rpud) {
393 ERROR_exit("NULL rpud!!!\n"
394 "rpud %p, ud %p\n", rpud, ud);
395 }
396 if( val == NULL ){
397 if (rpud->verb > 1) {
398 INFO_message("First call npts=%d\n", npts);
399 Show_RP_UD(rpud, "Top of init:\n");
400 }
401 if (rpud->pmeth != FFT_PHASE) {
402 ERROR_message("This should not happen");
403 }
404 if( npts > 0 ){ /* the "start notification" */
405 ncall = 0 ;
406 st = Dir2Type(rpud->dir);
407 if (st < 0) {
408 ERROR_message("Bad rpud->dir, assuming POLAR stimulus");
409 st = POL;
410 }
411 if (rpud->Dsamp) {
412 WARNING_message(
413 "Does not deal with slice timing offset correctly. Needs fixing");
414 }
415 /* nfft ? */
416 if (rpud->nfft == 0) rpud->nfft = csfft_nextup(rpud->nvals);
417 if (rpud->verb > 1) {
418 INFO_message(
419 "Data length = %d ; FFT length = %d; st %d, direc %d (%s)",
420 rpud->nvals,rpud->nfft,
421 st, rpud->dir, Phase_Dirs_lbl(rpud->dir)) ;
422 }
423
424 if( rpud->ftap > 0.0f ) {
425 xtap = mri_setup_taper( rpud->nvals , rpud->ftap ) ;
426 }
427
428 /* freq. step */
429 rpud->fstep = 1.0f/(rpud->nfft*rpud->dt);
430
431 /* response frequency */
432 rpud->Fresp[st] =
433 rpud->Fstim[st]*(float)rpud->n[st];
434
435 if (rpud->Fresp[st] <= 0.0f) {
436 ERROR_message("Bad rpud->Fresp !");
437 return;
438 }
439
440 /* allocate for fft array */
441 comp_array = (complex *) calloc( sizeof(complex) , rpud->nfft);
442 mag = (float *) calloc( sizeof(float) , rpud->nfft);
443 phz = (float *) calloc( sizeof(float) , rpud->nfft);
444 nzfreq = (byte *) calloc( sizeof(byte) , rpud->nfft);
445 stimharm = (int *) calloc( sizeof(int) , rpud->nfft);
446
447 /* signal bin */
448 if (!SetFreqBin(rpud->Fresp[st],
449 rpud->fstep, rpud->stk, rpud->stw,
450 rpud->nfft)) {
451 ERROR_message("Failed to set bins. Fresp=%f, "
452 "fstep=%f, stk=[%d,%d], stw=[%f,%f], nfft=%d\n"
453 "Is your dataset's TR of %f sec valid?",
454 rpud->Fresp[st], rpud->fstep,
455 rpud->stk[0], rpud->stk[1],
456 rpud->stw[0], rpud->stw[1],rpud->nfft,
457 rpud->dt);
458 return;
459 }
460 {
461 int nharm = 1, stk[2];
462 float fharm=0.0, stw[2];
463 while (( (fharm = nharm*rpud->Fresp[st]) -
464 rpud->nfft/2.0*rpud->fstep ) < -0.00001) {
465 /* Difference was added to get around round off errors
466 which might get you a frequency at the precipice as
467 in the bug reported by Phil Burton Nov. 2014 */
468 if (!SetFreqBin(fharm, rpud->fstep, stk, stw, rpud->nfft)) {
469 ERROR_message("Failed to set bins. Fresp=%f, fdiff=%f "
470 "fstep=%f, stk=[%d,%d], stw=[%f,%f], nfft=%d\n"
471 "Is your dataset's TR of %f sec valid?",
472 rpud->Fresp[st], fharm-rpud->nfft/2.0*rpud->fstep,
473 rpud->fstep,
474 rpud->stk[0], rpud->stk[1],
475 rpud->stw[0], rpud->stw[1],rpud->nfft,rpud->dt);
476 break;
477 }
478 stimharm[stk[0]] = nharm;
479 stimharm[stk[1]] = nharm;
480 if (rpud->verb > 2)
481 INFO_message("Freq. indices [%d(%.3f) %d(%.3f)] is harm %d"
482 " (fharm=%f. lim=%f)\n",
483 stk[0], stw[0], stk[1], stw[1], nharm,
484 fharm, rpud->nfft/2.0*rpud->fstep);
485 ++nharm;
486 }
487 }
488 /* mark frequencies to be used for noise variance estimate */
489 N_nzfreq=0;
490 for( jj=0 ; jj < rpud->nfft/2 ; jj++ ) {
491 if ( jj > rpud->stk[1] /* above freq. */
492 && !stimharm[jj] ) { /* not harmonic */
493 nzfreq[jj] = 1; ++N_nzfreq;
494 if (rpud->verb > 2) INFO_message("Freq. index %d is noise\n", jj);
495 }
496 }
497 rpud->dof[0] = 2; /* assuming estimate of phase and amp at one freq.*/
498 rpud->dof[1] = 2*N_nzfreq; /* phase and amp from freq. of noise */
499
500 /* init output sets */
501 if (rpud->spectra) {
502 for (ii=0; ii<2; ++ii) {
503 THD_3dim_dataset *oset=NULL;
504 char stmp[10+strlen(rpud->prefix)];
505
506 if (ii) {
507 if (rpud->verb > 1) INFO_message("Init. amp");
508 sprintf(stmp,"%s.%s.amp%s",
509 rpud->prefix, Phase_Dirs_lbl(rpud->dir),rpud->oext);
510 } else {
511 if (rpud->verb > 1) INFO_message("Init. phz");
512 sprintf(stmp,"%s.%s.phz%s",
513 rpud->prefix, Phase_Dirs_lbl(rpud->dir),rpud->oext);
514 }
515 oset = EDIT_empty_copy( rpud->iset ) ;
516 EDIT_dset_items( oset ,
517 ADN_prefix , stmp,
518 ADN_datum_all, MRI_float ,
519 ADN_nvals , rpud->nfft ,
520 ADN_ntt , rpud->nfft ,
521 ADN_none ) ;
522
523 DSET_UNMSEC(rpud->iset) ;
524 if( DSET_TIMEUNITS(rpud->iset) == UNITS_SEC_TYPE ){
525 EDIT_dset_items( oset ,
526 ADN_tunits , UNITS_HZ_TYPE ,
527 ADN_ttdel , rpud->fstep ,
528 ADN_nsl , 0 ,
529 ADN_none ) ;
530 } else {
531 WARNING_message("Units not seconds?");
532 }
533 for( jj=0 ; jj < rpud->nfft ; jj++ ) {
534 EDIT_substitute_brick( oset , jj , MRI_float , NULL ) ;
535 }
536 if (ii) {
537 rpud->amp = oset; oset = NULL;
538 } else {
539 rpud->phz = oset; oset = NULL;
540 }
541 }
542 } else {
543 rpud->amp = rpud->phz = NULL;
544 }
545
546 if (rpud->verb >1) {
547 Show_RP_UD(rpud, "@ End of init:\n");
548 }
549 } else { /* the "end notification" */
550 if (rpud->verb > 1) {
551 INFO_message("Last call\n");
552 }
553 st = -1;
554 if (xtap) free(xtap); xtap = NULL;
555 if (comp_array) free(comp_array); comp_array = NULL;
556 if (mag) free(mag); mag = NULL;
557 if (phz) free(phz); mag = NULL;
558 if (nzfreq) free(nzfreq); nzfreq = NULL;
559 }
560 return ;
561 }
562 if (rpud->verb > 3) {
563 INFO_message("call %d\n", ncall);
564 }
565 /* Now do the FFT */
566 for( jj=0 ; jj < rpud->nvals ; jj++ ) {
567 comp_array[jj].r = ts[jj]; comp_array[jj].i = 0.0f ;
568 }
569
570 if( xtap != NULL ){
571 for( jj=0 ; jj < rpud->nvals ; jj++ ){
572 comp_array[jj].r *= xtap[jj] ; comp_array[jj].i *= xtap[jj] ;
573 }
574 }
575 for( jj=rpud->nvals ; jj < rpud->nfft ; jj++ )
576 comp_array[jj].r = comp_array[jj].i = 0.0f ; /* zero pad */
577
578 csfft_cox( -1 , rpud->nfft, comp_array ) ; /* DFT */
579
580 /* Calculate phase, magnitude, etc. */
581 nzpow = 0.0;
582 preoff = rpud->pre * rpud->Fresp[st] * 2.0 * PI;
583 /* pre stim offset */
584
585 for( jj=0 ; jj < rpud->nfft ; jj++ ) {
586 mag[jj] = CABS(comp_array[jj]) ;
587 if (nzfreq[jj]) {
588 nzpow += (mag[jj]*mag[jj]);
589 }
590 phz[jj] = atan2(comp_array[jj].i, comp_array[jj].r);
591
592 /* remove offset due to pre stimulus period */
593 if (preoff >= 0) {
594 phz[jj] -= preoff;
595 if (phz[jj] < -PI) phz[jj] = 2.0*PI+phz[jj];
596 }
597 /* go from -pi ... pi to 0 ... 360 */
598 phz[jj] = PHASE_R2D(phz[jj]);
599 phz[jj] = PHASE_360(phz[jj]);
600 /* and change sign to reflect direction of stimulus */
601 if (rpud->dir < 0) phz[jj] = 360.0 - phz[jj];
602 /* Now divide by n so that phase is expressed in units of visual field */
603 if (rpud->n[st] > 1)
604 phz[jj] /= (float)rpud->n[st];
605
606 /* At this stage, phz[jj] is
607 theta+ if rpud->dir is +ve,
608 and theta- otherwise
609 (see labbook NIH-5, pp 103) */
610 }
611
612
613 /* Store the output vals */
614 /* take phase from the closest frequency with highest amplitude */
615 if (mag[rpud->stk[0]] > mag[rpud->stk[1]]){
616 val[0] = phz[rpud->stk[0]];
617 } else {
618 val[0] = phz[rpud->stk[1]];
619 }
620 /* linear interpolation between closest frequecies for amplitude */
621 val[1] = ( mag[rpud->stk[0]]*rpud->stw[0] +
622 mag[rpud->stk[1]]*rpud->stw[1] );
623 val[1] *= val[1]; /* square for power*/
624 val[1] /= (nzpow/(double)N_nzfreq); /* normalize by avg power of noise */
625 if (rpud->vox) {
626 jj=rpud->vox[ncall];
627 } else {
628 jj=ncall;
629 }
630 if (rpud->spectra) {
631 THD_insert_series(jj, rpud->amp, rpud->nfft, MRI_float, mag, 1);
632 THD_insert_series(jj, rpud->phz, rpud->nfft, MRI_float, phz, 1);
633 }
634 ncall++ ; return ;
635 }
636
DEL_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)637 static void DEL_tsfunc( double tzero, double tdelta ,
638 int npts, float ts[],
639 double ts_mean , double ts_slope ,
640 void *ud, int nbriks, float *val )
641 {
642 RP_UD *rpud = (RP_UD *)ud;
643 static int ncall = 0, reverse=0;
644 static float scale=0.0;
645 char * label=NULL; /* string containing stat. summary of results */
646 int errcode=0, ixyz=0;
647 float slp=0.0, delu=0.0, del=0.0, xcor=0.0, xcorCoef=0.0,vts=0.0,
648 vrvec=0.0, dtx=0.0;
649 static int st = -1;
650
651 if (!rpud) {
652 ERROR_exit("NULL rpud!!!\n"
653 "rpud %p, ud %p\n", rpud, ud);
654 }
655
656 /* Now intialize */
657 if( val == NULL ){
658 if (rpud->verb > 1) {
659 INFO_message("First call npts=%d\n", npts);
660 Show_RP_UD(rpud, "Top of init:\n");
661 }
662 if (rpud->pmeth != HILB_DELAY) {
663 ERROR_message("This should not happen");
664 }
665
666 if( npts > 0 ){ /* the "start notification" */
667 ncall = 0 ;
668 st = Dir2Type(rpud->dir);
669 if (st < 0) {
670 ERROR_message("Bad rpud->dir, assuming POLAR stimulus");
671 st = POL;
672 }
673
674 if (rpud->nfft != 0) { /* should allow for this and taper, in future */
675 ERROR_message("Control of nfft not allowed for hilbert delay");
676 return;
677 }
678 if (rpud->verb > 1) {
679 INFO_message("Data length = %d (npts=%d); st %d, direc %d (%s)",
680 rpud->nvals, npts,
681 st, rpud->dir, Phase_Dirs_lbl(rpud->dir)) ;
682 }
683 /* response frequency */
684 rpud->Fresp[st] =
685 rpud->Fstim[st]*(float)rpud->n[st];
686
687 if (rpud->Fresp[st] <= 0.0f) {
688 ERROR_message("Bad rpud->Fresp !");
689 return;
690 }
691 if (!rpud->rvec || rpud->rvec_num < 1) {
692 ERROR_message("No reference time series");
693 return;
694 }
695 if (rpud->rvec_len != rpud->nvals) {
696 ERROR_message( "Reference time series has %d vals, "
697 "time series has %d", rpud->rvec_len, rpud->nvals);
698 return;
699 }
700 if (rpud->Dsamp) {
701 WARNING_message(
702 "Does not deal with slice timing offset correctly. Needs fixing");
703 }
704 rpud->dof[0] = 2; /* two fit params*/
705 rpud->dof[1] = rpud->ort_adj; /* number of orts */
706 if (rpud->dof[1] < 2) rpud->dof[1] = 2; /* always linear detrend...*/
707
708 if (rpud->dir == CLW || rpud->dir == EXP) reverse = 1;
709 else reverse = 0;
710
711 /* initialize hilbert */
712 set_delay_verb(rpud->verb);
713 hilbertdelay_V2reset();
714
715 if (rpud->verb > 1) {
716 Show_RP_UD(rpud, "@ End of init:\n");
717 }
718 } else { /* the "end notification" */
719 if (rpud->verb > 1) {
720 INFO_message("Last call\n");
721 }
722 st = -1;
723 }
724 return ;
725 }
726
727 if (rpud->verb > 3) {
728 INFO_message("call %d\n", ncall);
729 }
730
731 /* get slice offset time
732 WARNING: NOT ALLOWING for offset from time series cropping*/
733 if (rpud->Dsamp) {
734 dtx = (float) (tzero / tdelta);
735 } else {
736 dtx = 0.0;
737 }
738
739 if (rpud->vox) {
740 ixyz=rpud->vox[ncall];
741 } else {
742 ixyz=ncall;
743 }
744
745 errcode =
746 hilbertdelay_V2 (ts, /* voxel time series */
747 rpud->rvec[rpud->iref], /* reference ts*/
748 npts, /* length of time series */
749 1, 0, /* Num. of segments and percent overlap */
750 1 ,0, /* not cleanup mode, no detrend */
751 dtx, /* timing offset */
752 1, /* remove bias */
753 &delu,&slp,&xcor,&xcorCoef,&vts,&vrvec);
754
755 if (errcode == 0) { /* If there are no errors, proceed */
756 hunwrap (delu, 1.0/rpud->dt, 1.0/rpud->Fresp[st], slp,
757 0, METH_DEGREES, reverse,
758 1.0/rpud->n[st], &del );
759
760 } else if (errcode == ERROR_LONGDELAY) {
761 if (0) {
762 WARNING_message("Errcode LONGDELAY at voxel %d\n", ixyz);
763 }
764
765 del = 0.0; /* Set all the variables to Null and don't set xcorCoef
766 to an impossible value*/
767 xcorCoef = 0.0;/* because the data might still be OK */
768 xcor = 0.0;
769
770 } else if (errcode == ERROR_QUIT) {
771 exit(1);
772 } else {
773 if (0) {
774 WARNING_message("Errcode %d at voxel %d\n", errcode, ixyz);
775 }
776
777 del = 0.0; /* Set all the variables to Null and set xcorCoef
778 to an impossible value*/
779 xcorCoef = NOWAYXCORCOEF;
780 xcor = 0.0;
781 }
782
783
784 /*----- Save results for this voxel -----*/
785 val[0] = del;
786 val[1] = xcorCoef;
787
788
789
790 ncall++ ; return ;
791 }
792
793
MaskSetup(THD_3dim_dataset * old_dset,THD_3dim_dataset * mask_dset,RP_UD * rpud,byte * cmask,int * ncmask,float mask_bot,float mask_top,int * mcount)794 byte *MaskSetup(THD_3dim_dataset *old_dset, THD_3dim_dataset *mask_dset,
795 RP_UD *rpud, byte *cmask, int *ncmask,
796 float mask_bot, float mask_top, int *mcount)
797 {
798 byte *mmm=NULL;
799 int ii=0, kk=0;
800
801 /* ------------- Mask business -----------------*/
802 if( mask_dset == NULL ){
803 mmm = NULL ;
804 if( rpud->verb )
805 INFO_message("%d voxels in the entire dataset (no mask)\n",
806 DSET_NVOX(old_dset)) ;
807 } else {
808 if( DSET_NVOX(mask_dset) != DSET_NVOX(old_dset) )
809 ERROR_exit("Input and mask datasets are not same dimensions!\n");
810 mmm = THD_makemask( mask_dset , 0 , mask_bot, mask_top ) ;
811 *mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
812 if( *mcount <= 0 ) {
813 ERROR_message("No voxels in the mask!\n") ;
814 return(NULL);
815 }
816 if( rpud->verb ) INFO_message("%d voxels in the mask dset %s\n",
817 *mcount, DSET_PREFIX(mask_dset)) ;
818 DSET_delete(mask_dset) ;
819 }
820
821 if( cmask != NULL ){
822 if( *ncmask != DSET_NVOX(old_dset) )
823 ERROR_exit("Input and cmask datasets are not same dimensions!\n");
824 if( mmm != NULL ){
825 for( ii=0 ; ii < DSET_NVOX(old_dset) ; ii++ )
826 mmm[ii] = (mmm[ii] && cmask[ii]) ;
827 free(cmask) ;
828 *mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
829 if( *mcount <= 0 ) {
830 ERROR_message("No voxels in the mask+cmask!\n") ;
831 return(NULL);
832 }
833 if( rpud->verb )
834 INFO_message("%d voxels in the mask+cmask\n",*mcount) ;
835 } else {
836 mmm = cmask ;
837 *mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
838 if( *mcount <= 0 ) {
839 ERROR_message("No voxels in the cmask!\n") ;
840 return(NULL);
841 }
842 if( rpud->verb ) INFO_message("%d voxels in the cmask\n",*mcount) ;
843 }
844 }
845
846 if (mmm) {
847 rpud->nmask=*mcount;
848 rpud->vox = (int *)calloc(sizeof(int),*mcount) ;
849 kk=0;
850 for (ii=0; ii < DSET_NVOX(old_dset) ; ii++ ) {
851 if (mmm[ii]) { rpud->vox[kk]=ii; ++kk; }
852 }
853 }
854
855 return(mmm);
856 }
857
main(int argc,char * argv[])858 int main( int argc , char * argv[] )
859 {
860 THD_3dim_dataset **new_dset=NULL, *old_dset=NULL,
861 *mask_dset=NULL, *idset=NULL;
862 char stmp[1024], **in_name=NULL;
863 int iarg=1 , ii, jj, kk, ll, nvox, nvals=1, isfloat=0, nset=0, stype=0;
864 int datum = MRI_float, mcount = 0, ncmask=0, dtused = 0;
865 byte *mmm=NULL, *cmask=NULL;
866 float mask_bot=1.0 , mask_top=-1.0,
867 *x0=NULL, *x1=NULL, *d0=NULL, *d1=NULL, *best=NULL ;
868 RP_UD rpud;
869
870 in_name = (char **)calloc(sizeof(char*),k_N);
871 new_dset = (THD_3dim_dataset **)calloc(sizeof(THD_3dim_dataset*),k_N);
872
873 rpud.pmeth = FFT_PHASE;
874 rpud.Dsamp = 0;
875 rpud.rvec = NULL; rpud.rvec_num = 0; rpud.rvec_len = 0; rpud.iref = 0;
876 rpud.ftap = 0.0f;
877 rpud.dt = 0.0f;
878 rpud.Fstim[ECC] = rpud.Fstim[POL] = 0.0f;
879 rpud.Fresp[ECC] = rpud.Fresp[POL] = 0.0f;
880 rpud.nfft = 0;
881 rpud.verb = 0;
882 rpud.iset = NULL;
883 rpud.nmask=0;
884 rpud.vox = NULL;
885 rpud.amp = NULL;
886 rpud.phz = NULL;
887 rpud.prefix = "RetinoPhase";
888 rpud.oext=NULL;
889 rpud.dir = NOT_SET;
890 rpud.n[ECC] = rpud.n[POL] = 1;
891 rpud.spectra = 0;
892 rpud.pre = 0.0f;
893 rpud.fixsum = 1;
894 rpud.ort_adj = 0;
895
896 /* rpud. = ; */
897
898
899 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
900 printf("Usage: 3dRetinoPhase [-prefix ppp] dataset\n"
901 " where dataset is a time series from a retinotpy stimulus\n"
902 "\n"
903 " -exp EXP: These four options specify the type of retinotpy \n"
904 " -con CON: stimulus. EXP and CON are for expanding and \n"
905 " -clw CLW : contracting rings, respectively. CLW and CCW are\n"
906 " -ccw CCW: for clockwise and counter clockwise moving polar\n"
907 " polar angle mapping stimuli. You can specify one, \n"
908 " or all stimuli in one command. When all are specified\n"
909 " polar angle stimuli, and eccentricity stimuli of \n"
910 " opposite directions are combined.\n"
911 " -prefix PREF: Prefix of output datasets. \n"
912 " PREF is suffixed with the following:\n"
913 " .ecc+ for positive (expanding) eccentricity (EXP)\n"
914 " .ecc- for negative (contracting) eccentricity (CON)\n"
915 " .pol+ for clockwise polar angle mapping (CLW)\n"
916 " .pol- for counterclockwise polar angle mapping (CCW)\n"
917 " At a minimum each input gets a phase dataset output. It contains\n"
918 " response phase (or delay) in degrees.\n"
919 " If both directions are given for polar and/or eccentricity\n"
920 " then a visual field angle data set is created.\n"
921 " The visual field angle is obtained by averaging phases of opposite\n"
922 " direction stimuli. The hemodynamic offset is half the phase difference.\n"
923 "\n"
924 " Each output also contains a thresholding sub-brick. Its type \n"
925 " depends on the phase estimation method (-phase_estimate).\n"
926 "\n"
927 " Note on the thresholding sub-bricks\n"
928 " -----------------------------------\n"
929 " Both FFT and DELAY values of -phase_estimate produce thresholding \n"
930 " sub-bricks with the phase estimates. Those thresholds have associated \n"
931 " significance levels, but they should be taken with a grain of \n"
932 " salt. There is no correction for autocorrelation, so the DOFs \n"
933 " are generous.\n"
934 " The program also attaches a thresholding sub-brick to the\n"
935 " visual field angle datasets which are estimated by averaging the phase\n"
936 " estimates in order to remove the hemodynamic offset. This composite \n"
937 " thresholding sub-brick contains at each voxel/node, the maximum\n"
938 " threshold from the datasets of stimli of opposite direction.\n"
939 " This thresholding sub-brick is for convenience, allowing you to\n"
940 " threshold with a mask that is the union of the individual\n"
941 " thresholded maps. Significance levels are purposefully not\n"
942 " attached. I don't know how to compute them properly.\n"
943 "\n"
944 " -spectra: Output amplitude and phase spectra datasets.\n"
945 " -Tstim T: Period of stimulus in seconds. This parameter does\n"
946 " not depend on the number of wedges or rings (Nr/Nw).\n"
947 " It is the duration of a full cycle of the stimulus.\n"
948 " Use -Tpol TPOL, and -Tecc TECC, to specify periods\n"
949 " for each stimulus type separately. -Tstim sets both \n"
950 " periods to T.\n"
951 " -nrings Nr: Nr is the number of rings in the stimulus. \n"
952 " The default is 1.\n"
953 " -nwedges Nw: Nw is the number of wedges in the stimulus. \n"
954 " The default is 1.\n"
955 " -ort_adjust: Number of DOF lost in detrending outside of this \n"
956 " program.\n"
957 " -pre_stim PRE: Blank period, in seconds, before stimulus began \n"
958 " -sum_adjust y/n: Adjust sum of angles for wrapping based on the\n"
959 " angle difference. Default is 'y'\n"
960 " -phase_estimate METH: Select method of phase estimation\n"
961 " METH == FFT uses the phase of the fundamental frequency.\n"
962 " METH == DELAY uses the 3ddelay approach for estimating\n"
963 " the phase. This requires the use of option\n"
964 " -ref_ts . See references [3] and [4] below. \n"
965 " The DELAY option appears to be good as the FFT for high SNR\n"
966 " and high duty cycle. See results produced by @Proc.PK.All_D\n"
967 " in the demo archive AfniRetinoDemo.tgz.\n"
968 " However,the DELAY option seems much better for low duty cycle stimuli.\n"
969 " It is not set as the default for backward compatibility. Positive and \n"
970 " negative feedback about this option are welcome.\n"
971 "\n"
972 " Thanks to Ikuko Mukai and Masaki Fukunaga for making the case \n"
973 " for DELAY's addition; they were right. \n"
974 "\n"
975 " -ref_ts REF_TS: 0 lag reference time series of response. This is\n"
976 " needed for the DELAY phase estimation method.\n"
977 " With the DELAY method, the phase results are comparable to \n"
978 " what you'd get with the following 3ddelay command:\n"
979 " For illustration, say you have stimuli of 32 second periods\n"
980 " with the polar stimuli having two wedges. After creating \n"
981 " the reference time series with waver (32 sec. block period \n"
982 " eccentricity, 32/2=16 sec. block period for polar), run \n"
983 " 4 3ddelay commands as such:\n"
984 " for an expanding ring of 32 second period:\n"
985 " 3ddelay -input exp.niml.dset \\\n"
986 " -ideal_file ECC.1D \\\n"
987 " -fs 0.5 -T 32 \\\n"
988 " -uD -nodsamp \\\n"
989 " -phzreverse -phzscale 1.0 \\\n"
990 " -prefix ecc+.del.niml.dset\\n"
991 "\n"
992 " Repeat for contracting ring, remove -phzreverse \n"
993 "\n"
994 " for clockwise two wedge of 32 second period:\n"
995 " 3ddelay -input clw.niml.dset \\\n"
996 " -ideal_file POL.1D \\\n"
997 " -fs 0.5 -T 16 \\\n"
998 " -uD -nodsamp \\\n"
999 " -phzreverse -phzscale 0.5 \\\n"
1000 " -prefix pol+.del.niml.dset\\n"
1001 "\n"
1002 " Repeat for counterclockwise remove -phzreverse \n"
1003 " Instead of the 3ddelay mess, all you do is run 3dRetinoPhase with the \n"
1004 " following extra options: "
1005 " -phase_estimate DELAY -ref_ts ECC.1D\n"
1006 " or -phase_estimate DELAY -ref_ts POL.1D\n"
1007 "\n"
1008 " If you are not familiar with the use of program 'waver' for creating\n"
1009 " reference time series, take a look at demo script @Proc.PK.All_D in\n"
1010 " AfniRetinoDemo.tgz.\n"
1011 "\n"
1012 " -multi_ref_ts MULTI_REF_TS: Multiple 0 lag reference time series. \n"
1013 " This allows you to test multiple regressors.\n"
1014 " The program will run a separate analysis for \n"
1015 " each regressor (column), and combine the results\n"
1016 " in the output dataset this way:\n"
1017 " ([.] denotes output sub-brick)\n"
1018 " [0]: Phase from regressor that yields the highest correlation coeff.\n"
1019 " [1]: Maximum correlation coefficient.\n"
1020 " [2]: Number of regressor that yields the highest correlation coeff.\n"
1021 " Counting begins at 1 (not 0)\n"
1022 " [3]: Phase from regressor 1\n"
1023 " [4]: Correlation coefficient from regressor 1\n"
1024 " [5]: Phase from regressor 2\n"
1025 " [6]: Correlation coefficient from regressor 2\n"
1026 " ... etc.\n"
1027 " In general, for regressor k (k starts at 1)\n"
1028 " [2*k+1] contains the Phase and [2*k+2] the Correlation coefficient\n"
1029 "\n"
1030 " N.B: If MULTI_REF_TS has only one timeseries, -multi_ref_ts produces\n"
1031 " an output identical to that of -ref_ts. \n"
1032 "\n"
1033 " See usage in @RetinoProc and demo data in\n"
1034 " https://afni.nimh.nih.gov/pub/dist/tgz/AfniRetinoDemo.tgz \n"
1035 "\n"
1036 "References for this program:\n"
1037 " [1] RW Cox. AFNI: Software for analysis and visualization of functional\n"
1038 " magnetic resonance neuroimages. \n"
1039 " Computers and Biomedical Research, 29: 162-173, 1996.\n"
1040 " [2] Saad Z.S., et al. SUMA: An Interface For Surface-Based Intra- And\n"
1041 " Inter-Subject Analysis With AFNI.\n"
1042 " Proc. 2004 IEEE International Symposium on Biomedical Imaging, 1510-1513\n"
1043 " If you use the DELAY method:\n"
1044 " [3] Saad, Z.S., et al. Analysis and use of FMRI response delays. \n"
1045 " Hum Brain Mapp, 2001. 13(2): p. 74-93.\n"
1046 " [4] Saad, Z.S., E.A. DeYoe, and K.M. Ropella, Estimation of FMRI \n"
1047 " Response Delays. Neuroimage, 2003. 18(2): p. 494-504.\n"
1048 "\n"
1049 /* options left over from 3dDFT.c , or not yet tested
1050 " -abs == output float dataset = abs(DFT)\n"
1051 " -nfft N == use 'N' for DFT length (must be >= #time points)\n"
1052 " -taper f == taper 'f' fraction of data at ends (0 <= f <= 1).\n"
1053 " [Hamming 'raised cosine' taper of f/2 of the ]\n"
1054 " [data length at each end; default is no taper]\n"
1055 " -slice_time_adjust y/n: Adjust for slice timing difference.\n"
1056
1057 Bring them back after testing. */
1058 ) ;
1059 PRINT_COMPILE_DATE ; exit(0) ;
1060 }
1061
1062 mainENTRY("3dRetinoPhase main"); machdep(); AFNI_logger("3dRetinoPhase",argc,argv);
1063 /* AUTHOR("ZIAD") ; */
1064 #ifdef USING_MCW_MALLOC
1065 enable_mcw_malloc() ;
1066 #endif
1067
1068 /*-- options --*/
1069 set_obliquity_report(0); /* silence obliquity */
1070
1071 while( iarg < argc && argv[iarg][0] == '-' ){
1072
1073 if( strncmp(argv[iarg],"-mask",5) == 0 ){
1074 if( mask_dset != NULL )
1075 ERROR_exit("Cannot have two -mask options!\n") ;
1076 if( iarg+1 >= argc )
1077 ERROR_exit("-mask option requires a following argument!\n");
1078 mask_dset = THD_open_dataset( argv[++iarg] ) ;
1079 if( mask_dset == NULL )
1080 ERROR_exit("Cannot open mask dataset!\n") ;
1081 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex )
1082 ERROR_exit("Cannot deal with complex-valued mask dataset!\n");
1083 iarg++ ; continue ;
1084 }
1085
1086 if( strncmp(argv[iarg],"-mrange",5) == 0 ){
1087 if( iarg+2 >= argc )
1088 ERROR_exit("-mrange option requires 2 following arguments!\n");
1089 mask_bot = strtod( argv[++iarg] , NULL ) ;
1090 mask_top = strtod( argv[++iarg] , NULL ) ;
1091 if( mask_top < mask_bot )
1092 ERROR_exit("-mrange inputs are illegal!\n") ;
1093 iarg++ ; continue ;
1094 }
1095
1096 if( strncmp(argv[iarg],"-verb",5) == 0 ){
1097 if( iarg+1 >= argc )
1098 ERROR_exit("-verb option requires 1 integer!\n");
1099 rpud.verb = (int)strtod( argv[++iarg] , NULL ) ;
1100 iarg++ ; continue ;
1101 }
1102
1103 if( strcmp(argv[iarg],"-cmask") == 0 ){ /* 16 Mar 2000 */
1104 if( iarg+1 >= argc )
1105 ERROR_exit("-cmask option requires a following argument!\n");
1106 cmask = EDT_calcmask( argv[++iarg] , &ncmask, 0 ) ;
1107 if( cmask == NULL ) ERROR_exit("Can't compute -cmask!\n");
1108 iarg++ ; continue ;
1109 }
1110
1111 if( strcmp(argv[iarg],"-taper") == 0 ){
1112 if( iarg+1 >= argc )
1113 ERROR_exit("-taper option requires an argument!\n");
1114 rpud.ftap = (float)strtod(argv[++iarg],NULL) ;
1115 if( rpud.ftap < 0.0f || rpud.ftap > 1.0f )
1116 ERROR_exit("Illegal value after -taper: %g",rpud.ftap) ;
1117 iarg++ ; continue ;
1118 }
1119
1120 if( strcmp(argv[iarg],"-sum_adjust") == 0 ){
1121 if( iarg+1 >= argc )
1122 ERROR_exit("-sum_adjust requires an argument\n");
1123 ++iarg;
1124 if (argv[iarg][0] == 'y' || argv[iarg][0] == 'Y') {
1125 rpud.fixsum = 1;
1126 } else if (argv[iarg][0] == 'n' || argv[iarg][0] == 'N') {
1127 rpud.fixsum = 0;
1128 } else {
1129 ERROR_exit("Illegal value (%s) after -sum_adjust\n", argv[iarg]);
1130 }
1131 iarg++ ; continue ;
1132 }
1133
1134 if( strcmp(argv[iarg],"-slice_time_adjust") == 0 ){
1135 if( iarg+1 >= argc )
1136 ERROR_exit("-slice_time_adjust requires an argument\n");
1137 ++iarg;
1138 if (argv[iarg][0] == 'y' || argv[iarg][0] == 'Y') {
1139 rpud.Dsamp = 1;
1140 } else if (argv[iarg][0] == 'n' || argv[iarg][0] == 'N') {
1141 rpud.Dsamp = 0;
1142 } else {
1143 ERROR_exit("Illegal value (%s) after -slice_time_adjust\n",
1144 argv[iarg]);
1145 }
1146 iarg++ ; continue ;
1147 }
1148
1149 if( strcmp(argv[iarg],"-phase_estimate") == 0 ){
1150 if( iarg+1 >= argc )
1151 ERROR_exit("-phase_estimate requires an argument\n");
1152 ++iarg;
1153 if (!strcmp(argv[iarg],"FFT") || !strcmp(argv[iarg],"fft")) {
1154 rpud.pmeth = FFT_PHASE;
1155 } else if ( !strcmp(argv[iarg],"3DDELAY") ||
1156 !strcmp(argv[iarg],"3ddelay") ||
1157 !strcmp(argv[iarg],"delay") ||
1158 !strcmp(argv[iarg],"DELAY")) {
1159 rpud.pmeth = HILB_DELAY;
1160 } else {
1161 ERROR_exit("Illegal value (%s) after -phase_estimate\n",
1162 argv[iarg]);
1163 }
1164 iarg++ ; continue ;
1165 }
1166
1167 if( strcmp(argv[iarg],"-ref_ts") == 0 ||
1168 strcmp(argv[iarg],"-multi_ref_ts") == 0){
1169 MRI_IMAGE * flim=NULL;
1170 float *fp=NULL;
1171 int iv=0;
1172 if( iarg+1 >= argc )
1173 ERROR_exit("-ref_ts/-multi_ref_ts require an argument\n");
1174 ++iarg;
1175 if (!(flim = mri_read_1D(argv[iarg]))) {
1176 ERROR_exit("Illegal value (%s) after -phase_estimate\n",
1177 argv[iarg]);
1178 }
1179 rpud.rvec_len = flim->nx;
1180 rpud.rvec_num = flim->ny;
1181 if (flim->ny != 1) {
1182 if (!strcmp(argv[iarg-1],"-ref_ts")) {
1183 ERROR_exit("Only one column can be present in -ref_ts file.\n"
1184 "Have %d columns in %s.\n"
1185 "If you want to run the multi_ref version use\n"
1186 "-multi_ref_ts option instead, but understand that\n"
1187 "the output datasets have a differing content.\n"
1188 ,flim->ny, argv[iarg]);
1189 } else {
1190 /* multiref mode */
1191 }
1192 }
1193 rpud.rvec = (float **)calloc(flim->ny, sizeof(float *));
1194 fp = MRI_FLOAT_PTR(flim);
1195 for (iv=0; iv<flim->ny; ++iv) {
1196 rpud.rvec[iv] = (float *)calloc(flim->nx, sizeof(float));
1197 for (ii=0; ii<flim->nx; ++ii) {
1198 rpud.rvec[iv][ii] = *fp; ++fp;
1199 }
1200 }
1201 mri_free(flim); flim = NULL;
1202 iarg++ ; continue ;
1203 }
1204
1205
1206 if( strcmp(argv[iarg],"-fstim") == 0 ){
1207 rpud.Fstim[ECC] = rpud.Fstim[POL] = (float)strtod(argv[++iarg],NULL) ;
1208 if( rpud.Fstim[ECC] <= 0.0f )
1209 ERROR_exit("Illegal frequency after -fstim: %g",rpud.Fstim[ECC]) ;
1210 iarg++ ; continue ;
1211 }
1212
1213 if( strcmp(argv[iarg],"-fecc") == 0 ){
1214 rpud.Fstim[ECC] = (float)strtod(argv[++iarg],NULL) ;
1215 if( rpud.Fstim[ECC] <= 0.0f )
1216 ERROR_exit("Illegal frequency after -fecc: %g",rpud.Fstim[ECC]) ;
1217 iarg++ ; continue ;
1218 }
1219 if( strcmp(argv[iarg],"-fpol") == 0 ){
1220 rpud.Fstim[POL] = (float)strtod(argv[++iarg],NULL) ;
1221 if( rpud.Fstim[POL] <= 0.0f )
1222 ERROR_exit("Illegal frequency after -fpol: %g",rpud.Fstim[POL]) ;
1223 iarg++ ; continue ;
1224 }
1225
1226 if( strcmp(argv[iarg],"-Tstim") == 0 ){
1227 rpud.Fstim[ECC]= rpud.Fstim[POL]= 1.0/(float)strtod(argv[++iarg],NULL) ;
1228 if( rpud.Fstim[ECC] <= 0.0f )
1229 ERROR_exit("Illegal period after -Tstim: %g",1.0/rpud.Fstim[ECC]) ;
1230 iarg++ ; continue ;
1231 }
1232
1233 if( strcmp(argv[iarg],"-Tecc") == 0 ){
1234 rpud.Fstim[ECC]= 1.0/(float)strtod(argv[++iarg],NULL) ;
1235 if( rpud.Fstim[ECC] <= 0.0f )
1236 ERROR_exit("Illegal period after -Tecc: %g",1.0/rpud.Fstim[ECC]) ;
1237 iarg++ ; continue ;
1238 }
1239
1240 if( strcmp(argv[iarg],"-Tpol") == 0 ){
1241 rpud.Fstim[POL]= 1.0/(float)strtod(argv[++iarg],NULL) ;
1242 if( rpud.Fstim[POL] <= 0.0f )
1243 ERROR_exit("Illegal period after -Tpol: %g",1.0/rpud.Fstim[POL]) ;
1244 iarg++ ; continue ;
1245 }
1246
1247 if( strncmp(argv[iarg],"-pre_stim", 5) == 0 ){
1248 if( iarg+1 >= argc )
1249 ERROR_exit("-pre_stim requires an argument\n");
1250
1251 rpud.pre = (float)strtod(argv[++iarg],NULL) ;
1252 if( rpud.pre < 0.0f )
1253 ERROR_exit("Illegal value after -pre_stim: %g",rpud.pre) ;
1254 iarg++ ; continue ;
1255 }
1256
1257 if( strncmp(argv[iarg],"-nw",3) == 0 ){
1258 if( iarg+1 >= argc )
1259 ERROR_exit("-nwedges requires an integer\n");
1260 rpud.n[POL] = (int)strtod(argv[++iarg],NULL) ;
1261 if( rpud.n[POL] <1 || rpud.n[POL] > 10 )
1262 ERROR_exit("Bad value of %d for -nwedges. "
1263 "Should be a small positive integer.",rpud.n[POL]) ;
1264 iarg++ ; continue ;
1265 }
1266 if( strncmp(argv[iarg],"-nr",3) == 0 ){
1267 if( iarg+1 >= argc )
1268 ERROR_exit("-nrings requires an integer\n");
1269 rpud.n[ECC] = (int)strtod(argv[++iarg],NULL) ;
1270 if( rpud.n[ECC] <1 || rpud.n[ECC] > 10 )
1271 ERROR_exit("Bad value of %d for -nrings. "
1272 "Should be a small positive integer.",rpud.n[ECC]) ;
1273 iarg++ ; continue ;
1274 }
1275
1276 if( strcmp(argv[iarg],"-dir") == 0 ){
1277 ERROR_exit("-dir now obsolete, use -exp, -con, -clw, or -ccw\n");
1278 if( iarg+1 >= argc )
1279 ERROR_exit("-dir option requires an argument!\n");
1280 ++iarg;
1281 if (!strncmp(argv[iarg],"exp", 3)) rpud.dir = EXP;
1282 else if (!strncmp(argv[iarg],"con", 3) ||
1283 !strncmp(argv[iarg],"cnt", 3)) rpud.dir = CONT;
1284 else if (!strncmp(argv[iarg],"cw", 2) ||
1285 !strncmp(argv[iarg],"clo", 3) ||
1286 !strncmp(argv[iarg],"clw", 3)) rpud.dir = CLW;
1287 else if (!strncmp(argv[iarg],"ccw", 3) ||
1288 !strncmp(argv[iarg],"cou", 3)) rpud.dir = CCW;
1289 else if (!strncmp(argv[iarg],"+", 1) ||
1290 !strncmp(argv[iarg],"pos", 3)) rpud.dir = POS;
1291 else if (!strncmp(argv[iarg],"-", 1) ||
1292 !strncmp(argv[iarg],"neg", 3)) rpud.dir = NEG;
1293 else {
1294 ERROR_exit("Illegal value %s after -dir\n", argv[iarg]);
1295 }
1296 iarg++ ; continue ;
1297 }
1298
1299 if( strcmp(argv[iarg],"-input") == 0 ){
1300 ERROR_exit("-input now obsolete, use -exp, -con, -cw, or -ccw\n");
1301
1302 if( iarg+1 >= argc )
1303 ERROR_exit("-input option requires a dset\n");
1304 in_name[0] = argv[++iarg] ;
1305 iarg++ ; continue ;
1306 }
1307
1308 if( strcmp(argv[iarg],"-exp") == 0 ){
1309 if( iarg+1 >= argc )
1310 ERROR_exit("-exp option requires a dset\n");
1311 in_name[k_EXP] = argv[++iarg] ;
1312 iarg++ ; continue ;
1313 }
1314
1315 if( strcmp(argv[iarg],"-con") == 0 ){
1316 if( iarg+1 >= argc )
1317 ERROR_exit("-con option requires a dset\n");
1318 in_name[k_CON] = argv[++iarg] ;
1319 iarg++ ; continue ;
1320 }
1321
1322 if( strcmp(argv[iarg],"-cw") == 0 ||
1323 strcmp(argv[iarg],"-clw") == 0){
1324 if( iarg+1 >= argc )
1325 ERROR_exit("-clw (or -cw) option requires a dset\n");
1326 in_name[k_CLW] = argv[++iarg] ;
1327 iarg++ ; continue ;
1328 }
1329
1330 if( strcmp(argv[iarg],"-ccw") == 0 ){
1331 if( iarg+1 >= argc )
1332 ERROR_exit("-ccw option requires a dset\n");
1333 in_name[k_CCW] = argv[++iarg] ;
1334 iarg++ ; continue ;
1335 }
1336
1337 if( strcmp(argv[iarg],"-prefix") == 0 ){
1338 rpud.prefix = argv[++iarg] ;
1339 if( !THD_filename_ok(rpud.prefix) )
1340 ERROR_exit("-prefix %s is illegal!",rpud.prefix) ;
1341 iarg++ ; continue ;
1342 }
1343
1344 if( strcmp(argv[iarg],"-ort_adjust") == 0 ){
1345 rpud.ort_adj += atoi(argv[++iarg]);
1346 iarg++ ; continue ;
1347 }
1348
1349 if( strcmp(argv[iarg],"-spectra") == 0 ){
1350 rpud.spectra = 1 ; iarg++ ; continue ;
1351 }
1352
1353 if( strcmp(argv[iarg],"-detrend") == 0 ){
1354 dtused = 1;
1355 iarg++ ; continue ;
1356 }
1357
1358 if( strcmp(argv[iarg],"-nfft") == 0 ){
1359 rpud.nfft = (int)strtod(argv[++iarg],NULL) ;
1360 if( rpud.nfft <= 2 ){
1361 WARNING_message("Illegal -nfft value on command line") ;
1362 rpud.nfft = 0 ;
1363 } else {
1364 ii = csfft_nextup(rpud.nfft) ;
1365 if( ii > rpud.nfft ){
1366 WARNING_message("Replacing -nfft=%d with next largest"
1367 "legal value=%d",
1368 rpud.nfft,ii) ;
1369 rpud.nfft = ii ;
1370 }
1371 }
1372 iarg++ ; continue ;
1373 }
1374
1375 ERROR_exit("ILLEGAL option: %s\n",argv[iarg]) ;
1376 }
1377
1378 if (dtused && rpud.verb) {
1379 INFO_message("Linear detrend is always on, -detrend is useless.") ;
1380 }
1381
1382 if( !in_name )
1383 ERROR_exit("No datasets on command line!?") ;
1384
1385 if (rpud.verb && rpud.rvec && rpud.rvec_num > 1) {
1386 INFO_message("Multi-ref mode with %d reference time series\n"
1387 , rpud.rvec_num);
1388 }
1389
1390 nset = 0;
1391 for (stype=0; stype<k_N; ++stype) {
1392 if (in_name[stype]) {
1393 old_dset = THD_open_dataset( in_name[stype] ) ;
1394 if( !ISVALID_DSET(old_dset) )
1395 ERROR_exit("Can't open dataset %s\n",in_name[stype]);
1396
1397
1398 if( DSET_NVALS(old_dset) < 2 )
1399 ERROR_exit("Can't use dataset with < 2 values per voxel!\n") ;
1400
1401 if( DSET_NUM_TIMES(old_dset) < 2 ){
1402 WARNING_message("Input dataset is not 3D+time; assuming TR=1.0") ;
1403 EDIT_dset_items( old_dset ,
1404 ADN_ntt , DSET_NVALS(old_dset) ,
1405 ADN_ttorg , 0.0 ,
1406 ADN_ttdel , 1.0 ,
1407 ADN_tunits , UNITS_SEC_TYPE ,
1408 NULL ) ;
1409 }
1410
1411 if (!nset) {
1412 rpud.nvals = DSET_NUM_TIMES(old_dset);
1413 if (rpud.nvals < 10) {
1414 ERROR_exit( "Dataset has too few (%d) time points.\n",
1415 rpud.nvals);
1416 }
1417 rpud.dt = DSET_TR_SEC(old_dset);
1418 if (rpud.dt > 100) {
1419 ERROR_exit("TR of %f sec of %s is very suspicously high.\n"
1420 "Program assumes this is a mistake in the header.\n"
1421 "You can use 3drefit's -TR option to fix the problem.\n",
1422 rpud.dt, in_name[stype]);
1423 }
1424 mmm = MaskSetup(old_dset, mask_dset, &rpud, cmask,
1425 &ncmask, mask_bot, mask_top, &mcount);
1426 } else {
1427 if (DSET_NUM_TIMES(old_dset) != rpud.nvals) {
1428 ERROR_exit( "Dataset %s has %d time points while %s has %d\n",
1429 in_name[stype], DSET_NUM_TIMES(old_dset),
1430 DSET_PREFIX(rpud.iset), rpud.nvals);
1431 }
1432 if (DSET_TR_SEC(old_dset) != rpud.dt) {
1433 ERROR_exit( "Dataset %s has a TR of %f sec while %s has %f sec\n",
1434 in_name[stype], DSET_TR_SEC(old_dset),
1435 DSET_PREFIX(rpud.iset), rpud.dt);
1436 }
1437 }
1438 rpud.iset = old_dset;
1439 rpud.dir = Phase_Type_to_Dir(stype);
1440
1441 /*------------- ready to compute new dataset -----------*/
1442 if (rpud.verb)
1443 INFO_message("Processing %s, %d time points\n",
1444 in_name[stype], DSET_NVALS(old_dset));
1445
1446
1447 rpud.oext = "";
1448 if (STRING_HAS_SUFFIX(rpud.prefix, ".niml.dset")) {
1449 rpud.oext = ".niml.dset";
1450 rpud.prefix[strlen(rpud.prefix)-strlen(rpud.oext)] = '\0';
1451 } else if (STRING_HAS_SUFFIX(rpud.prefix, ".1D")) {
1452 rpud.oext = ".1D";
1453 rpud.prefix[strlen(rpud.prefix)-strlen(rpud.oext)] = '\0';
1454 } else if (STRING_HAS_SUFFIX(rpud.prefix, ".1D.dset")) {
1455 rpud.oext = ".1D.dset";
1456 rpud.prefix[strlen(rpud.prefix)-strlen(rpud.oext)] = '\0';
1457 } else if (STRING_HAS_SUFFIX(in_name[stype], ".niml.dset")) {
1458 rpud.oext = ".niml.dset";
1459 } else if (STRING_HAS_SUFFIX(in_name[stype], ".1D") ||
1460 STRING_HAS_SUFFIX(in_name[stype], ".1D.dset")) {
1461 rpud.oext = ".1D.dset";
1462 }
1463 sprintf(stmp,"%s.%s%s",
1464 rpud.prefix, Phase_Dirs_lbl(rpud.dir),rpud.oext);
1465
1466 if (rpud.verb > 1) INFO_message("Output time\n");
1467 if (rpud.pmeth == FFT_PHASE) {
1468 new_dset[stype] = MAKER_4D_to_typed_fbuc(
1469 old_dset , /* input dataset */
1470 stmp , /* output prefix */
1471 datum , /* output datum */
1472 0 , /* ignore count */
1473 1 , /* linear detrend */
1474 2 , /* number of briks */
1475 RP_tsfunc, /* timeseries processor */
1476 (void *)(&rpud), /* data for tsfunc */
1477 mmm,
1478 0 /* Allow auto scaling of output */
1479 ) ;
1480
1481 if( new_dset[stype] != NULL ){
1482 sprintf(stmp,"Phz@%.3fHz", rpud.Fresp[Dir2Type(rpud.dir)]);
1483 EDIT_BRICK_LABEL(new_dset[stype], 0, stmp);
1484 sprintf(stmp,"PwR@%.3fHz", rpud.Fresp[Dir2Type(rpud.dir)]);
1485 EDIT_BRICK_LABEL(new_dset[stype], 1, stmp);
1486 EDIT_BRICK_TO_FIFT(new_dset[stype],1,rpud.dof[0],rpud.dof[1]);
1487 } else {
1488 ERROR_exit("Unable to compute output dataset!\n") ;
1489 }
1490 } else if (rpud.pmeth == HILB_DELAY){
1491 for (rpud.iref = 0; rpud.iref < rpud.rvec_num; ++rpud.iref) {
1492 if (!(idset = MAKER_4D_to_typed_fbuc(
1493 old_dset , /* input dataset */
1494 stmp , /* output prefix */
1495 datum , /* output datum */
1496 0 , /* ignore count */
1497 1 , /* linear detrend */
1498 2 , /* number of briks */
1499 DEL_tsfunc, /* timeseries processor */
1500 (void *)(&rpud), /* data for tsfunc */
1501 mmm,
1502 0 /* Allow auto scaling of output */
1503 ))) {
1504 ERROR_exit("Unable to compute %d output dataset!\n",
1505 rpud.iref) ;
1506 }
1507 if (rpud.verb) {
1508 INFO_message("Done with delays for %d/%d\n",
1509 rpud.iref, rpud.rvec_num);
1510 }
1511 if (!new_dset[stype]) {
1512 /* first, pass, take output as is */
1513 new_dset[stype] = idset;
1514 /* fix labels */
1515 if (rpud.rvec_num == 1) {
1516 sprintf(stmp,"Phz_Delay");
1517 EDIT_BRICK_LABEL(new_dset[stype], 0, stmp);
1518 sprintf(stmp,"Corr.Coef.");
1519 EDIT_BRICK_LABEL(new_dset[stype], 1, stmp);
1520 EDIT_BRICK_TO_FICO(new_dset[stype], 1, rpud.rvec_len,
1521 rpud.dof[0],rpud.dof[1]);
1522 } else {
1523 if (rpud.verb > 1) {
1524 INFO_message("Adding Summary results\n", rpud.rvec_num);
1525 }
1526 sprintf(stmp,"Phz_Delay@Max.Corr.");
1527 EDIT_BRICK_LABEL(new_dset[stype], 0, stmp);
1528 sprintf(stmp,"Max.Corr.Coef.");
1529 EDIT_BRICK_LABEL(new_dset[stype], 1, stmp);
1530 /* add the best phase winner */
1531 best = (float *)malloc(DSET_NVOX(idset)*sizeof(float));
1532 for (ii=0; ii<DSET_NVOX(idset); ++ii) best[ii] = 1;
1533 EDIT_add_brick(new_dset[stype], MRI_float,
1534 0.0, best);
1535 EDIT_BRICK_LABEL(new_dset[stype], 2, "best_fit_regressor");
1536 /* the bulk of idset is now in the summary output
1537 make a copy of it because we'll also append it below
1538 Note that only idset's bricks are being transferred,
1539 so there is some leakage each time a new idset
1540 gets generated without wiping the old one clean...*/
1541 idset = EDIT_full_copy(idset,"gingrich");
1542 }
1543 }
1544 if (rpud.rvec_num > 1) { /* have multiple series */
1545 if (rpud.verb > 1) {
1546 INFO_message("Adding %d/%d results\n",
1547 rpud.iref+1,rpud.rvec_num);
1548 }
1549 /* add phase */
1550 EDIT_add_brick(new_dset[stype], DSET_BRICK_TYPE(idset,0), 0.0,
1551 DSET_ARRAY(idset,0));
1552 sprintf(stmp,"Phz_Delay@%02d", rpud.iref+1);
1553 EDIT_BRICK_LABEL(new_dset[stype],
1554 DSET_NVALS(new_dset[stype])-1, stmp);
1555 /* add threshold */
1556 EDIT_add_brick(new_dset[stype], DSET_BRICK_TYPE(idset,1), 0.0,
1557 DSET_ARRAY(idset,1));
1558 sprintf(stmp,"Corr.Coef.@%02d", rpud.iref+1);
1559 EDIT_BRICK_LABEL(new_dset[stype],
1560 DSET_NVALS(new_dset[stype])-1, stmp);
1561 EDIT_BRICK_TO_FICO(new_dset[stype],
1562 DSET_NVALS(new_dset[stype])-1,
1563 rpud.rvec_len,
1564 rpud.dof[0],rpud.dof[1]);
1565
1566 /* swap values in 1st two sets with the max. */
1567 if (rpud.iref > 0) {
1568 d0 = (float *)DSET_ARRAY(new_dset[stype], 0);
1569 x0 = (float *)DSET_ARRAY(new_dset[stype], 1);
1570 best = (float *)DSET_ARRAY(new_dset[stype], 2);
1571 d1 = (float *)DSET_ARRAY(new_dset[stype],
1572 DSET_NVALS(new_dset[stype])-2);
1573 x1 = (float *)DSET_ARRAY(new_dset[stype],
1574 DSET_NVALS(new_dset[stype])-1);
1575 for (ii=0; ii<DSET_NVOX(idset); ++ii) {
1576 if (x1[ii] > x0[ii]) {
1577 x0[ii] = x1[ii];
1578 d0[ii] = d1[ii];
1579 best[ii] = rpud.iref+1;
1580 }
1581 }
1582 }
1583 }
1584 }
1585 }
1586
1587
1588 /* write the output */
1589 if (rpud.verb) {
1590 INFO_message("Writing delay results");
1591 }
1592 tross_Copy_History( old_dset , new_dset[stype] ) ;
1593 tross_Make_History("3dRetinoPhase", argc, argv ,new_dset[stype]) ;
1594 /* refresh stats, what with the constant updating */
1595 DSET_KILL_STATS(new_dset[stype]); THD_load_statistics(new_dset[stype]) ;
1596 DSET_write( new_dset[stype] ) ;
1597 WROTE_DSET( new_dset[stype] ) ;
1598
1599 if (rpud.amp) {
1600 tross_Copy_History( old_dset , rpud.amp ) ;
1601 tross_Make_History( "3dRetinoPhase" , argc, argv , rpud.amp ) ;
1602 for (ii=0; ii<rpud.nfft; ++ii) {
1603 sprintf(stmp,"Amp@%.3fHz", ii*rpud.fstep);
1604 EDIT_BRICK_LABEL(rpud.amp, ii, stmp);
1605 }
1606 DSET_write( rpud.amp ) ;
1607 WROTE_DSET( rpud.amp ) ;
1608 DSET_delete(rpud.amp); rpud.amp = NULL;
1609 } else if (rpud.spectra) {
1610 ERROR_exit("Unable to compute output dataset!\n") ;
1611 }
1612 if (rpud.phz) {
1613 tross_Copy_History( old_dset , rpud.phz ) ;
1614 tross_Make_History( "3dRetinoPhase" , argc, argv , rpud.phz ) ;
1615 for (ii=0; ii<rpud.nfft; ++ii) {
1616 sprintf(stmp,"Phz@%.3fHz", ii*rpud.fstep);
1617 EDIT_BRICK_LABEL(rpud.phz, ii, stmp);
1618 }
1619 DSET_write( rpud.phz ) ;
1620 WROTE_DSET( rpud.phz ) ;
1621 DSET_delete(rpud.phz); rpud.phz = NULL;
1622 } else if (rpud.spectra) {
1623 ERROR_exit("Unable to compute output dataset!\n") ;
1624 }
1625
1626 DSET_delete(old_dset); old_dset = NULL;
1627 ++nset;
1628 } /* if stimulus type given */
1629 } /* for each stimulus type */
1630
1631 /* Now see if you can do combos */
1632 if (new_dset[k_EXP] && new_dset[k_CON]) {
1633 THD_3dim_dataset *cset=NULL;
1634 if (rpud.verb) {
1635 INFO_message("Combining bidirectional eccentricity data");
1636 }
1637 rpud.dir = CONT; /* not really, but sign does not matter at this point */
1638 cset = Combine_Opposites(new_dset[k_EXP], new_dset[k_CON],
1639 &rpud);
1640 tross_Copy_History( cset , new_dset[k_EXP] ) ;
1641 DSET_write( cset ) ;
1642 WROTE_DSET( cset ) ;
1643 DSET_delete(cset); cset = NULL;
1644 }
1645
1646 if (new_dset[k_CLW] && new_dset[k_CCW]) {
1647 THD_3dim_dataset *cset=NULL;
1648 if (rpud.verb) {
1649 INFO_message("Combining bidirectional polar data");
1650 }
1651 rpud.dir = CLW; /* not really, but sign does not matter at this point */
1652 cset = Combine_Opposites(new_dset[k_CLW], new_dset[k_CCW],
1653 &rpud);
1654 tross_Copy_History( cset , new_dset[k_CLW] ) ;
1655 DSET_write( cset ) ;
1656 WROTE_DSET( cset ) ;
1657 DSET_delete(cset); cset = NULL;
1658 }
1659
1660 if (rpud.verb>1) {
1661 INFO_message("Cleanup");
1662 }
1663 for (stype=0; stype<k_N; ++stype)
1664 if (new_dset[stype]) DSET_delete(new_dset[stype]);
1665
1666 free(in_name);
1667 free(new_dset);
1668
1669 exit(0) ;
1670 }
1671