1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "afni.h"
8 #include "afni_plugin.h"
9
10 #ifndef ALLOW_PLUGINS
11 # error "Plugins not properly set up -- see machdep.h"
12 #endif
13
14 /* -----------------------START---------------------------------*/
15 /* definition and decleration part to suport the main algorithm */
16
17 #include <stdlib.h>
18 #include <stdio.h>
19 #include <math.h>
20
21
22 /*------------------ macros to return workspace at exit -------------------*/
23
24 #define ZFREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
25
26 #undef ZFREE_WORKSPACE
27 #define ZFREE_WORKSPACE \
28 do{ ZFREEUP(bptr) ; ZFREEUP(sptr) ; ZFREEUP(fptr) ; \
29 ZFREEUP(cptr) ; ZFREEUP(fxar) ; ZFREEUP(fac) ; \
30 ZFREEUP(dtr) ; \
31 } while(0)
32 /*-------------------------------------------------------------------------*/
33
34
35 /* Definitions of prototypes and declaration of support functions
36 this is taken from the list of include files that I use in the original code*/
37
38
39 /*-------------------------------------------------------------------*/
40 /* COMPLEX STRUCTURE */
41
42 typedef struct {
43 float x, y, z;
44 } fXYZ;
45
46
47 /***********************************************************************
48 Plugin to extract 3D+time time courses whos index or xyz corrdinates
49 match a certain criterion
50 ************************************************************************/
51 typedef struct
52 {
53 int nxx; /* number of voxels in the x direction */
54 int nyy; /* number of voxels in the y direction */
55 int nzz; /* number of voxels in the z direction */
56 char *dsetname; /* prefix of data set analyzed */
57 int ignore; /* number ofpoints to ignore from time courses */
58 int ln; /* length of FMRI vector */
59 int dtrnd; /* remove linear trend or just the mean */
60 int errcode; /* error code number returned from function */
61 int out; /* flag for writing to a file */
62 int outts; /* flag for writing time series to a file */
63 int format;
64 int iloc;
65 int xloc;
66 int yloc;
67 int zloc;
68 int ncols;
69 int nrows;
70 float * indvect; /* vector that will hold the list of indices */
71 fXYZ * xyzvect; /* vecor that will hold the list of xyz coordinates */
72 char * strout;
73 char * strin;
74 FILE * outwritets;
75 FILE * outlogfile;
76 char outname[PLUGIN_MAX_STRING_RANGE]; /* output data file name */
77 }extract_data;
78
79 /*--------------------- string to 'help' the user --------------------*/
80
81 static char helpstring[] =
82 " 4D Dump Plugin\n\n"
83 "This plugin allows the extraction of selected FMRI time series into an ascii file. \n"
84 "The time series to be extracted are specified either by their AFNI index or \n"
85 "by their AFNI XYZ corrdinates.\n\n"
86 "Plugin Inputs\n"
87 " 1- Data :\n"
88 " 3D+time -> Chooser for 3D+time data set.\n"
89 " Ignore -> Specify the number of points to ignore from the beginning \n"
90 " of each voxel time series.\n"
91 " Dtrnd -> (Y/N) Apply detrending to time series. \n"
92 " If the dtrnd option is off, time series do not have zero mean.\n\n"
93 " 2- Mask file :\n"
94 " Mask File -> The index or X Y Z coordinates that specify wich voxel time series\n"
95 " should be output are specified in an ascii file.\n"
96 " Each line in the ascii file holds the index or X Y Z coordinates \n"
97 " of the voxel time series to be extracted. Each line of the ascii \n"
98 " file can hold many values, the user chooses which ones correspond \n"
99 " to i or the X Y Z triplet. The X Y Z coordinates used for masking\n"
100 " correspond to those shown in the AFNI window and NOT the graph window. \n"
101 " You can see these coordinates in the upper left side of the AFNI window.\n"
102 " To do so, you must switch the voxel coordinate units from mm to slice coordinates.\n"
103 " Define Datamode -> Misc -> Voxel Coords ?\n"
104 " Masking files could be the ascii output files of other plugins such as\n"
105 " '3Ddump' and 'Hilber Delay98'.\n"
106 " N Columns -> Total number of values on each line in the ascii mask file.\n\n"
107 " 3- Index Mask : (optional)\n"
108 " i col. -> Column number where the AFNI indices for the time series to be extracted \n"
109 " are stored in the ascii mask file.\n"
110 " 4- XYZ Mask : (optional)\n"
111 " x,y,z col. -> Specify the column numbers where X Y Z triplets are stored\n\n"
112 " 5- Output :\n"
113 " Filename -> Name of the ascii file where the time series are written.\n"
114 " If no output filename is specified, the ouput is the prefix\n"
115 " of the 3D+time input brick with the extenstion '.4Ddump' appended to it.\n"
116 " A LOG file, 'Filename.log' is also written to disk.\n"
117 " The log file contains all the parameters settings used\n"
118 " for generating 'Filename'.\n"
119 " Format -> ts[1] ts[2] .... -> one time series per line.\n"
120 " i x y z ts[1] ts[2] .... -> index, X Y Z coordinates and time series per line.\n\n"
121 " PS: If all the voxels specified in the mask file exist in the data file then each line\n"
122 " in the output file correspond to the line in the mask file. Otherwise you should use\n"
123 " the option i x y z ts[1] ts[2] .... to know which time series came from which voxel\n\n"
124 "If you have/find questions/comments/bugs about the plugin, \n"
125 "send me an E-mail: ziad@image.bien.mu.edu\n\n"
126 " Ziad Saad Jan 6 97, lastest update Feb 23 98.\n\n"
127 ;
128
129 /*--------------------- strings for output format --------------------*/
130
131 static char * yn_strings[] = { "n" , "y" };
132 static char * format_strings[] = { "i x y z ts[1] ..." , "ts[1] ts[2] ..." };
133
134 #define NUM_YN_STRINGS (sizeof(yn_strings)/sizeof(char *))
135 #define NUM_FORMAT_STRINGS (sizeof(yn_strings)/sizeof(char *))
136
137
138 #define YUP 1
139 #define NOPE 0
140
141 #define ERROR_NOSUCHVOXEL 1 /* Nothing to do in function */
142 #define ERROR_FILEREAD 2
143 #define ERROR_FILEWRITE 2
144 #define ERROR_OPTIONS 17
145 #define ERROR_OUTCONFLICT 19
146
147 /*----------------- prototypes for internal routines -----------------*/
148
149 static char * EXTRACT_main( PLUGIN_interface * ) ; /* the entry point */
150
151 static void EXTRACT_tsfunc( double T0 , double TR ,
152 int npts , float ts[] , double ts_mean , double ts_slope ,
153 void * udp) ;
154
155 static void show_ud (extract_data* ud);
156
157 static void write_ud (extract_data* ud);
158
159 static void indexTOxyz (extract_data* ud, int ncall, int *xpos , int *ypos , int *zpos);
160
161 static void error_report (extract_data* ud, int ncall );
162
163 static void writets (extract_data* ud,float * ts, int ncall);
164
165 static float * extract_index (char *fname, int ind_col_loc, int ncols, int *nrows, int *Err);
166
167 static fXYZ * extract_xyz (char *fname, int x_col_loc, int y_col_loc, int z_col_loc, int ncols, int *nrows, int *Err);
168
169 static void disp_vect (float *v,int l);
170
171 static int filexists (char *f_name);
172
173 static int f_file_size (char *f_name);
174
175 static int * PLUTO_4D_to_nothing (THD_3dim_dataset * old_dset , int ignore , int detrend ,
176 generic_func * user_func , void * user_data );
177
178 /*---------------------------- global data ---------------------------*/
179
180 static PLUGIN_interface * global_plint = NULL ;
181
182 /***********************************************************************
183 Set up the interface to the user:
184 1) Create a new interface using "PLUTO_new_interface";
185
186 2) For each line of inputs, create the line with "PLUTO_add_option"
187 (this line of inputs can be optional or mandatory);
188
189 3) For each item on the line, create the item with
190 "PLUTO_add_dataset" for a dataset chooser,
191 "PLUTO_add_string" for a string chooser,
192 "PLUTO_add_number" for a number chooser.
193 ************************************************************************/
194
195 DEFINE_PLUGIN_PROTOTYPE
196
PLUGIN_init(int ncall)197 PLUGIN_interface * PLUGIN_init( int ncall )
198 {
199 PLUGIN_interface * plint ; /* will be the output of this routine */
200
201 if( ncall > 0 ) return NULL ; /* only one interface */
202 CHECK_IF_ALLOWED("4DDUMP","4D Dump") ; /* 30 Sep 2016 */
203
204 /*---------------- set titles and call point ----------------*/
205
206 plint = PLUTO_new_interface( "4D Dump" ,
207 "Extract voxel time courses given their index or XYZ coordinates" ,
208 helpstring ,
209 PLUGIN_CALL_VIA_MENU , EXTRACT_main ) ;
210
211 global_plint = plint ; /* make global copy */
212
213 /*--------- 1st line: Input dataset and mask files ---------*/
214
215 PLUTO_add_option( plint ,
216 "Data" , /* label at left of input line */
217 "Data" , /* tag to return to plugin */
218 TRUE /* is this mandatory? */
219 ) ;
220
221 PLUTO_add_dataset( plint ,
222 "3D+time" , /* label next to button */
223 ANAT_ALL_MASK , /* take only EPI datasets */
224 FUNC_ALL_MASK , /* No fim funcs */
225 DIMEN_4D_MASK | /* need 3D+time datasets */
226 BRICK_ALLREAL_MASK /* need real-valued datasets */
227 ) ;
228
229 PLUTO_add_number( plint ,
230 "Ignore" , /* label next to chooser */
231 0 , /* smallest possible value */
232 50 , /* largest possible value (inactivated for now)*/
233 0 , /* decimal shift (none in this case) */
234 0 , /* default value */
235 FALSE /* allow user to edit value? */
236 ) ;
237 PLUTO_add_string( plint ,
238 "Dtrnd" , /* label next to textfield */
239 2,yn_strings, /* strings to choose among */
240 1 /* Default option */
241 ) ;
242
243 /*---------- 2nd line: Mask file info ----------*/
244 PLUTO_add_option( plint ,
245 "Mask file" , /* label at left of input line */
246 "Mask" , /* tag to return to plugin */
247 TRUE /* is this mandatory? */
248 ) ;
249
250 PLUTO_add_string( plint , "Mask File" , 0 , NULL , 19 ) ;
251
252 PLUTO_add_number( plint ,
253 "N Columns" , /* label next to chooser */
254 1 , /* smallest possible value */
255 1000 , /* largest possible value (inactivated for now)*/
256 0 , /* decimal shift (none in this case) */
257 3 , /* default value */
258 TRUE /* allow user to edit value? */
259 ) ;
260
261
262 /*---------- 3rd line: index mask location ----------*/
263
264 PLUTO_add_option( plint ,
265 "Index Mask ?" , /* label at left of input line */
266 "Index" , /* tag to return to plugin */
267 FALSE /* is this mandatory? */
268 ) ;
269
270 PLUTO_add_number( plint ,
271 "i col." , /* label next to chooser */
272 1 , /* smallest possible value */
273 1000 , /* largest possible value (inactivated for now)*/
274 0 , /* decimal shift (none in this case) */
275 1 , /* default value */
276 TRUE /* allow user to edit value? */
277 ) ;
278
279
280 /*---------- 4th line: xyz mask location ----------*/
281
282 PLUTO_add_option( plint ,
283 "XYZ Mask ?" , /* label at left of input line */
284 "XYZ" , /* tag to return to plugin */
285 FALSE /* is this mandatory? */
286 ) ;
287
288 PLUTO_add_number( plint ,
289 "x col." , /* label next to chooser */
290 1 , /* smallest possible value */
291 1000 , /* largest possible value */
292 0 , /* decimal shift (none in this case) */
293 2 , /* default value */
294 TRUE /* allow user to edit value? */
295 ) ;
296
297 PLUTO_add_number( plint ,
298 "y col." , /* label next to chooser */
299 1 , /* smallest possible value */
300 1000 , /* largest possible value */
301 0 , /* decimal shift (none in this case) */
302 3 , /* default value */
303 TRUE /* allow user to edit value? */
304 ) ;
305
306
307 PLUTO_add_number( plint ,
308 "z col." , /* label next to chooser */
309 1 , /* smallest possible value */
310 1000 , /* largest possible value */
311 0 , /* decimal shift (none in this case) */
312 4 , /* default value */
313 TRUE /* allow user to edit value? */
314 ) ;
315
316 /*---------- 5th line: output stuff ----------*/
317
318 PLUTO_add_option( plint ,
319 "Output" , /* label at left of input line */
320 "Output" , /* tag to return to plugin */
321 TRUE /* is this mandatory? */
322 ) ;
323
324
325 PLUTO_add_string( plint , "Filename" , 0 , NULL , 19 ) ;
326
327 PLUTO_add_string( plint ,
328 "Format" , /* label next to textfield */
329 2,format_strings, /* strings to choose among */
330 0 /* Default option */
331 ) ;
332
333 return plint ;
334 }
335
336 /***************************************************************************
337 Main routine for this plugin (will be called from AFNI).
338 If the return string is not NULL, some error transpired, and
339 AFNI will popup the return string in a message box.
340 ****************************************************************************/
341
EXTRACT_main(PLUGIN_interface * plint)342 static char * EXTRACT_main( PLUGIN_interface * plint )
343 {
344 extract_data uda,*ud;
345 MRI_IMAGE * tsim;
346 MCW_idcode * idc ; /* input dataset idcode */
347 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
348 char *tmpstr , * str ;
349 int ntime, nvec , Err=0 , itmp, nprfx;
350 float * vec , fs , T ;
351 char * tag; /* plugin option tag */
352
353 /* Allocate as much character space as Bob specifies in afni.h + a bit more */
354
355 tmpstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
356
357 if (tmpstr == NULL)
358 return "********************\n"
359 "Could not Allocate\n"
360 "a teeni weeni bit of\n"
361 "Memory ! \n"
362 "********************\n";
363
364 ud = &uda; /* ud now points to an allocated space */
365 ud->errcode = 0; /*reset error flag */
366
367 /*--------------------------------------------------------------------*/
368 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
369
370 /*--------- go to first input line ---------*/
371
372 tag = PLUTO_get_optiontag(plint) ;
373
374 if (tag == NULL)
375 {
376 return "************************\n"
377 "Bad 1st line option \n"
378 "************************" ;
379 }
380
381 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
382 old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
383 if( old_dset == NULL )
384 return "*************************\n"
385 "Cannot find Input Dataset\n"
386 "*************************" ;
387
388 ud->dsetname = DSET_PREFIX (old_dset);
389
390 ud->ignore = PLUTO_get_number(plint) ; /* get number item */
391
392 str = PLUTO_get_string(plint) ;
393 ud->dtrnd = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
394
395
396
397 /*--------- loop over remaining options ---------*/
398
399
400 ud->iloc = -1;
401 ud->xloc = -1;
402 ud->yloc = -1;
403 ud->zloc = -1;
404 do
405 {
406 tag = PLUTO_get_optiontag(plint) ;
407 if (tag == NULL) break;
408 if (strcmp (tag, "Mask") == 0)
409 {
410 ud->strin = PLUTO_get_string(plint) ;
411 ud->ncols = PLUTO_get_number(plint) ;
412 continue;
413 }
414
415 if (strcmp (tag, "Index") == 0)
416 {
417 ud->iloc = PLUTO_get_number(plint) ; /* get number item */
418 continue;
419 }
420
421 if (strcmp (tag, "XYZ") == 0)
422 {
423 ud->xloc = PLUTO_get_number(plint) ; /* get number item */
424 ud->yloc = PLUTO_get_number(plint) ; /* get number item */
425 ud->zloc = PLUTO_get_number(plint) ; /* get number item */
426 continue;
427 }
428
429 if (strcmp (tag, "Output") == 0)
430 {
431 ud->strout = PLUTO_get_string(plint) ;
432
433 str = PLUTO_get_string(plint) ;
434 ud->format = (int)PLUTO_string_index( str , NUM_FORMAT_STRINGS , format_strings );
435 continue;
436 }
437
438 } while (1);
439 /* ------------------ Check for some errorsor inconsistencies ------------- */
440
441 if (ud->iloc == -1 && ud->xloc == -1)
442 {
443 return "**************************\n"
444 "At least iloc or x/y/zloc\n"
445 "must be specified\n"
446 "**************************\n"
447 ;
448 }
449
450 if (ud->iloc != -1 && ud->xloc != -1)
451 {
452 return "***************************\n"
453 "iloc AND x/y/zloc can not\n"
454 "be simultaneously specified\n"
455 "***************************\n"
456 ;
457 }
458
459
460 /* ------------------Done with user parameters ---------------------------- */
461
462 /* Now loadup that index list or the xyz list */
463 if (ud->iloc != -1)
464 {
465 itmp = 0; /* might want to give option of setting it to number of rows if*/
466 /* the users know it, otherwise, it is automatically determined*/
467 ud->indvect = extract_index (ud->strin, ud->iloc, ud->ncols, &itmp, &Err);
468 }
469 else /* assuming the only other case is that x y z are specified */
470 {
471 itmp = 0;
472 ud->xyzvect = extract_xyz (ud->strin , ud->xloc , ud->yloc , ud->zloc , ud->ncols, &itmp, &Err);
473 }
474
475
476 ud->nrows = itmp;
477
478 switch (Err)
479 {
480 case (0):
481 break;
482 case (1):
483 return "****************************************\n"
484 "index location should be > 0 and < ncols\n"
485 "****************************************\n";
486 case (2):
487 return "********************************************\n"
488 "file size and number of columns do not match\n"
489 "********************************************\n";
490 case (3):
491 return "**********************\n"
492 "Can't find matrix file\n"
493 "**********************\n";
494 case (4):
495 return "*****************\n"
496 "ncols must be > 0\n"
497 "*****************\n";
498 case (5):
499 return "****************************************\n"
500 "x/y/z column numbers can NOT be the same\n"
501 "****************************************\n";
502 default:
503 return "****************************************\n"
504 "Should not have gotten here .....\n"
505 "****************************************\n";
506
507 }
508
509 /* check to see if the field is empty */
510 if (ud->strout == NULL) nprfx = 0;
511 else nprfx = 1;
512
513 /* check if the size is larger than 0. I did not want to check for this unless it's allocated */
514 if (nprfx == 1 && (int)strlen (ud->strout) == 0) nprfx = 0;
515
516
517 if (nprfx == 0) /* no output file is specified */
518 {
519 sprintf ( tmpstr , "%s.4Ddump" , DSET_PREFIX (old_dset));
520 ud->strout = tmpstr;
521 }
522
523
524 if (filexists(ud->strout) == 1)
525 {
526 return "*******************************\n"
527 "Outfile exists, won't overwrite\n"
528 "*******************************\n";
529 }
530 ud->outwritets = fopen (ud->strout ,"w");
531
532 sprintf ( tmpstr , "%s.log" , ud->strout);
533 if (filexists(tmpstr) == 1)
534 {
535 return "*******************************\n"
536 "Outfile.log exists, won't overwrite\n"
537 "*******************************\n";
538 }
539
540 ud->outlogfile = fopen (tmpstr ,"w");
541
542 if ((ud->outwritets == NULL) || (ud->outlogfile == NULL) )
543 {
544 ud->errcode = ERROR_FILEWRITE;
545
546 return "***********************\n"
547 "Could Not Write Outfile\n"
548 "***********************\n";
549 }
550
551
552 ud->nxx = (int)old_dset->daxes->nxx; /* get data set dimensions */
553 ud->nyy = (int)old_dset->daxes->nyy;
554 ud->nzz = (int)old_dset->daxes->nzz;
555
556 /* ready to dump the log file */
557 write_ud (ud);
558
559 /*------------- ready to compute new dataset -----------*/
560
561 PLUTO_4D_to_nothing (old_dset , ud->ignore , ud->dtrnd ,
562 EXTRACT_tsfunc , (void *)ud );
563
564 fclose (ud->outlogfile);
565 fclose (ud->outwritets);
566
567 if (tmpstr) free (tmpstr);
568 return NULL ; /* null string returned means all was OK */
569 }
570
571 /**********************************************************************
572 Function that does the real work
573 ***********************************************************************/
574
EXTRACT_tsfunc(double T0,double TR,int npts,float ts[],double ts_mean,double ts_slope,void * udp)575 static void EXTRACT_tsfunc( double T0 , double TR ,
576 int npts , float ts[] , double ts_mean , double ts_slope ,
577 void * udp)
578 {
579 static int nvox = -1, ncall = -1;
580 extract_data uda,*ud;
581 float xcor=0.0 , tmp=0.0 , tmp2 = 0.0 , dtx = 0.0 ,\
582 slp = 0.0 , vts = 0.0 , vrvec = 0.0 , rxcorCoef = 0.0;
583 int i , is_ts_null , status , opt , actv , zpos , ypos , xpos , hit;
584
585 ud = &uda;
586 ud = (extract_data *) udp;
587
588
589 /** is this a "notification"? **/
590
591
592 if( ncall == -1 || ncall == nvox)
593 {
594 if( npts > 0 ){ /* the "start notification" */
595 PLUTO_popup_meter( global_plint ) ; /* progress meter */
596 nvox = npts ; /* keep track of */
597 ncall = 0 ; /* number of calls */
598
599 } else { /* the "end notification" */
600 PLUTO_set_meter( global_plint , 100 ) ; /* set meter to 100% */
601
602 }
603 return ;
604 }
605
606
607 hit = 0;
608 ud->ln = npts;
609
610 if (ud->iloc != -1) /* for each ncall, find out if this index is wanted*/
611 {
612 for (i=0;i<ud->nrows;++i)
613 {
614 if (ud->indvect[i] == (float)ncall)
615 {
616 hit = 1;
617 writets (ud,ts,ncall);
618 i = ud->nrows;
619 }
620 }
621 }
622 else
623 {
624 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos);
625 for (i=0;i<ud->nrows;++i)
626 {
627 if (ud->xyzvect[i].x == (float)xpos)
628 {
629 if (ud->xyzvect[i].y == (float)ypos)
630 {
631 if (ud->xyzvect[i].z == (float)zpos)
632 {
633 hit = 1;
634 writets (ud,ts,ncall);
635 i = ud->nrows;
636 }
637 }
638 }
639 }
640
641 }
642
643
644 if (ud->errcode == 0) /* if there are no errors, proceed */
645 {/* */
646 }/* ud->errcode == 0 outer loop */
647
648 /** set the progress meter to the % of completion **/
649 ncall++ ;
650
651 PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
652
653 ud->errcode = 0; /* Reset error to nothing */
654
655 return ;
656 }
657
658 /* ************************************************************ */
659 /* function to display user data input (debugging stuff) */
660 /* ************************************************************ */
661
show_ud(extract_data * ud)662 static void show_ud (extract_data* ud)
663 {
664 printf ("\n\nUser Data Values at location :\n");
665 printf ("ud->dsetname= %s\n",ud->dsetname);
666 printf ("ud->strin= %s\n",ud->strin);
667 printf ("ud->strout= %s\n",ud->strout);
668 printf ("ud->nxx= %d\n",ud->nxx);
669 printf ("ud->nyy= %d\n",ud->nyy);
670 printf ("ud->nzz= %d\n",ud->nzz);
671 printf ("ud->iloc= %d\n",ud->iloc);
672 printf ("ud->xloc= %d\n",ud->xloc);
673 printf ("ud->yloc= %d\n",ud->yloc);
674 printf ("ud->zloc= %d\n",ud->zloc);
675 printf ("ud->ncols= %d\n",ud->ncols);
676 printf ("ud->nrows= %d\n",ud->nrows);
677 printf ("ud->ignore= %d\n",ud->ignore);
678 printf ("ud->dtrnd= %d\n",ud->dtrnd);
679 printf ("ud->format= %d\n",ud->format);
680 printf ("Hit enter to continue\a\n\n");
681 getchar ();
682 return;
683 }
684
685 /* ************************************************************ */
686 /* function to write user data input to log file */
687 /* ************************************************************ */
688
write_ud(extract_data * ud)689 static void write_ud (extract_data* ud)
690 {
691 fprintf (ud->outlogfile,"\n\nUser Data Values \n");
692 fprintf (ud->outlogfile,"ud->dsetname= %s\n",ud->dsetname);
693 fprintf (ud->outlogfile,"ud->strin= %s\n",ud->strin);
694 fprintf (ud->outlogfile,"ud->strout= %s\n",ud->strout);
695 fprintf (ud->outlogfile,"ud->nxx= %d\n",ud->nxx);
696 fprintf (ud->outlogfile,"ud->nyy= %d\n",ud->nyy);
697 fprintf (ud->outlogfile,"ud->nzz= %d\n",ud->nzz);
698 fprintf (ud->outlogfile,"ud->iloc= %d\n",ud->iloc);
699 fprintf (ud->outlogfile,"ud->xloc= %d\n",ud->xloc);
700 fprintf (ud->outlogfile,"ud->yloc= %d\n",ud->yloc);
701 fprintf (ud->outlogfile,"ud->zloc= %d\n",ud->zloc);
702 fprintf (ud->outlogfile,"ud->ncols= %d\n",ud->ncols);
703 fprintf (ud->outlogfile,"ud->nrows= %d\n",ud->nrows);
704 fprintf (ud->outlogfile,"ud->ignore= %d\n",ud->ignore);
705 fprintf (ud->outlogfile,"ud->dtrnd= %d\n",ud->dtrnd);
706 fprintf (ud->outlogfile,"ud->format= %d\n",ud->format);
707 fprintf (ud->outlogfile,"\nThe format for the output file is the following:\n");
708 switch (ud->format)
709 {
710 case (0):
711 fprintf (ud->outlogfile,"ncall\txpos\typos\tzpos\tts[0]\tts[1]\t...\n");
712 break;
713 case (1):
714 fprintf (ud->outlogfile,"ts[0]\tts[1]\t...\n");
715 break;
716
717 }
718
719 return;
720 }
721
722 /* ************************************************************ */
723 /* function to compute x, y, z coordinates from the index */
724 /* ************************************************************ */
725
indexTOxyz(extract_data * ud,int ncall,int * xpos,int * ypos,int * zpos)726 static void indexTOxyz (extract_data* ud, int ncall, int *xpos , int *ypos , int *zpos)
727 {
728 *zpos = (int)ncall / (int)(ud->nxx*ud->nyy);
729 *ypos = (int)(ncall - *zpos * ud->nxx * ud->nyy) / (int)ud->nxx;
730 *xpos = ncall - ( *ypos * ud->nxx ) - ( *zpos * ud->nxx * ud->nyy ) ;
731 return;
732 }
733
734 /* ************************************************************ */
735 /* function to report errors encountered to the logfile */
736 /* Only errors that happen during runtime are handeled here, */
737 /* the others are handeled instantaneously, and need not be */
738 /* logged */
739 /* ************************************************************ */
740
error_report(extract_data * ud,int ncall)741 void error_report (extract_data* ud, int ncall )
742 {
743 int xpos,ypos,zpos;
744
745 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos);
746
747 switch (ud->errcode)
748 {
749
750 default:
751 fprintf (ud->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",ud->errcode);
752 break;
753 }
754 fprintf (ud->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos );
755 return;
756 }
757
758 /* *************************************************************** */
759 /* function to write the time course into a line in the given file */
760 /* *************************************************************** */
761
writets(extract_data * ud,float * ts,int ncall)762 void writets (extract_data * ud,float * ts, int ncall)
763
764 {
765 int i,xpos,ypos,zpos;
766
767 switch (ud->format)
768 {
769 case (0):
770 indexTOxyz (ud,ncall,&xpos , &ypos , &zpos);
771 fprintf (ud->outwritets, "%d\t%d\t%d\t%d\t",ncall,xpos, ypos,zpos);
772 break;
773 case (1):
774 break;
775 default:
776 break;
777 }
778
779 for (i=0;i<ud->ln;++i)
780 {
781 fprintf (ud->outwritets, "%f\t",ts[i]);
782 }
783 fprintf (ud->outwritets,"\n");
784 }
785
786 /* *************************************************************** */
787 /* function to extract x y z colums form a matrix format file */
788 /* *************************************************************** */
extract_xyz(char * fname,int x_col_loc,int y_col_loc,int z_col_loc,int ncols,int * nrows,int * Err)789 static fXYZ * extract_xyz (char *fname, int x_col_loc, int y_col_loc, int z_col_loc, int ncols, int *nrows, int *Err)
790 {/*extract_xyz*/
791
792 float tmp, *tmpX;
793 int sz,i,indx,tst;
794 div_t divstuff,tempx,tempy,tempz;
795 FILE * INFILE;
796 fXYZ * xyzvect=NULL ;
797
798 /* ncols must be > 0 */
799 if (ncols <= 0)
800
801 {
802 *Err = 4;
803 return (xyzvect); /* ncols <= 0 !*/
804 }
805
806
807 /* ind_col_loc must be > 0 (1st column is 1, and it must be less than ncols) */
808 if (x_col_loc <= 0 || x_col_loc > ncols || y_col_loc <= 0 || y_col_loc > ncols || z_col_loc <= 0 || z_col_loc > ncols )
809 {
810 *Err = 1;
811 return (xyzvect); /* ?_col_loc should be >0 and <ncols*/
812 }
813
814 /* if the number of rows is given, compute required size */
815 if (*nrows > 0)
816 {
817 sz = *nrows * ncols;
818 }
819 else
820 { /* dtermine the size and compute the number of rows, check if it's an integer*/
821 sz = f_file_size (fname);
822 if (sz == -1)
823 {
824 *Err = 3;
825 return (xyzvect);
826 }
827 divstuff = div (sz,ncols);
828 if (divstuff.rem != 0)
829 {
830 *Err = 2;
831 return (xyzvect); /* size of file and number of columns don't match */
832 }
833 else
834 {
835 *nrows = divstuff.quot;
836 }
837 }
838
839 tst = (x_col_loc - y_col_loc) * (x_col_loc - z_col_loc) * (y_col_loc - z_col_loc);
840 if (tst == 0)
841 {
842 *Err = 5;
843 return (xyzvect); /* columns specified are the same */
844 }
845
846 /* Allocate and check for necessary space */
847 xyzvect = (fXYZ *) calloc (sz+2,sizeof(fXYZ));
848
849 if (xyzvect == NULL)
850 {
851 printf ("\nFatal Error : Failed to Allocate memory\a\n");
852 printf ("Abandon Lab Immediately !\n\n");
853 return NULL ;
854 };
855
856 INFILE = fopen (fname,"r");
857 if (INFILE == NULL)
858 {
859 *Err = 3; /* can't find file */
860 return (xyzvect);
861 }
862
863 tempx = div (x_col_loc,ncols);
864 tempy = div (y_col_loc,ncols);
865 tempz = div (z_col_loc,ncols);
866
867 for (i=0;i<sz;++i)
868 {
869 fscanf (INFILE,"%f",&tmp);
870 divstuff = div ((i+1),ncols);
871 if (divstuff.rem != 0)
872 indx = divstuff.quot;
873 else
874 indx = divstuff.quot - 1;
875
876 if (divstuff.rem == tempx.rem)
877 xyzvect[indx].x = tmp;
878
879 else if (divstuff.rem == tempy.rem)
880 xyzvect[indx].y = tmp;
881
882 else if (divstuff.rem == tempz.rem)
883 xyzvect[indx].z = tmp;
884 }
885
886
887 *Err = 0;
888 return (xyzvect);
889
890 }/*extract_xyz*/
891
892
893 /* *************************************************************** */
894 /* function to extract an indices columnform a matrix format file */
895 /* *************************************************************** */
896
897
extract_index(char * fname,int ind_col_loc,int ncols,int * nrows,int * Err)898 static float * extract_index (char *fname, int ind_col_loc, int ncols, int *nrows, int *Err)
899 {/*extract_index*/
900
901 float tmp, *indxvect=NULL;
902 int sz,i;
903 div_t divstuff,temp;
904 FILE * INFILE;
905
906 /* ncols must be > 0 */
907 if (ncols <= 0)
908
909 {
910 *Err = 4;
911 return (indxvect); /* ncols <= 0 !*/
912 }
913
914
915 /* ind_col_loc must be > 0 (1st column is 1, and it must be less than ncols) */
916 if (ind_col_loc <= 0 || ind_col_loc > ncols)
917 {
918 *Err = 1;
919 return (indxvect); /* ind_col_loc should be >0 and <ncols*/
920 }
921
922 /* if the number of rows is given, compute required size */
923 if (*nrows > 0)
924 {
925 sz = *nrows * ncols;
926 }
927 else
928 { /* dtermine the size and compute the number of rows, check if it's an integer*/
929 sz = f_file_size (fname);
930 if (sz == -1)
931 {
932 *Err = 3;
933 return (indxvect);
934 }
935 divstuff = div (sz,ncols);
936 if (divstuff.rem != 0)
937 {
938 *Err = 2;
939 return (indxvect); /* size of file and number of columns don't match */
940 }
941 else
942 {
943 *nrows = divstuff.quot;
944 }
945 }
946
947 /* Allocate and check for necessary space */
948 indxvect = (float *) calloc (sz+2,sizeof(float));
949
950 if (indxvect == NULL)
951 {
952 printf ("\nFatal Error : Failed to Allocate memory\a\n");
953 printf ("Abandon Lab Immediately !\n\n");
954 return NULL ;
955 };
956
957 INFILE = fopen (fname,"r");
958 if (INFILE == NULL)
959 {
960 *Err = 3; /* can't find file */
961 return (indxvect);
962 }
963
964 temp = div (ind_col_loc,ncols);
965
966 for (i=0;i<sz;++i)
967 {
968 fscanf (INFILE,"%f",&tmp);
969 divstuff = div ((i+1),ncols);
970 if (divstuff.rem == temp.rem)
971 {
972 if (divstuff.rem != 0) indxvect[divstuff.quot] = tmp;
973 else indxvect[divstuff.quot-1] = tmp;
974 }
975 }
976
977
978 *Err = 0;
979 return (indxvect);
980
981 }/*extract_index*/
982
983 /* *************************************************************** */
984 /* function to return the size of a float numbers file */
985 /* *************************************************************** */
986
f_file_size(char * f_name)987 static int f_file_size (char *f_name)
988
989 {
990
991
992 int cnt=0,ex;
993 float buf;
994
995 FILE*internal_file;
996
997 internal_file = fopen (f_name,"r");
998 if (internal_file == NULL) {
999 return (-1);
1000 }
1001 ex = fscanf (internal_file,"%f",&buf);
1002 while (ex != EOF)
1003 {
1004 ++cnt;
1005 ex = fscanf (internal_file,"%f",&buf);
1006 }
1007
1008
1009 fclose (internal_file);
1010 return (cnt);
1011 }
1012
1013 /* *************************************************************** */
1014 /* function to test if a file exists */
1015 /* *************************************************************** */
1016
filexists(char * f_name)1017 static int filexists (char *f_name)
1018 {/*filexists*/
1019 FILE *outfile;
1020
1021 outfile = fopen (f_name,"r");
1022 if (outfile == NULL)
1023 return (0);
1024 else
1025 fclose (outfile);
1026 return (1);
1027
1028 }/*filexists*/
1029
1030
1031 /* *************************************************************** */
1032 /* function to display a vector (debugging) */
1033 /* *************************************************************** */
1034
disp_vect(float * v,int l)1035 static void disp_vect (float *v,int l)
1036 {
1037 int i;
1038
1039 printf ("\n");
1040 if ((l-1) == 0)
1041 {
1042 printf ("V = %f\n",*v);
1043 }
1044 else
1045 {
1046 for (i=0;i<l;++i)
1047 {
1048 printf ("V[%d] = %f\t",i,v[i]);
1049 }
1050 printf ("\n");
1051 }
1052 return;
1053
1054 }
1055
1056 /* *************************************************************** */
1057 /* function to retrieve 4D data from AFNI */
1058 /* (adapted from BDW's function PLUTO_4D_to_typed_fith */
1059 /* *************************************************************** */
1060
PLUTO_4D_to_nothing(THD_3dim_dataset * old_dset,int ignore,int detrend,generic_func * user_func,void * user_data)1061 static int * PLUTO_4D_to_nothing (THD_3dim_dataset * old_dset , int ignore , int detrend ,
1062 generic_func * user_func, void * user_data )
1063 {
1064
1065 byte ** bptr = NULL ; /* one of these will be the array of */
1066 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
1067 float ** fptr = NULL ; /* (depending on input datum type) */
1068 complex ** cptr = NULL ;
1069
1070 float * fxar = NULL ; /* array loaded from input dataset */
1071 float * fac = NULL ; /* array of brick scaling factors */
1072 float * dtr = NULL ; /* will be array of detrending coeff */
1073
1074 float val , d0fac , d1fac , x0,x1;
1075 double tzero=0.0 , tdelta , ts_mean , ts_slope ;
1076 int ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox ;
1077 static int retval;
1078 register int kk ;
1079
1080 /*----------------------------------------------------------*/
1081 /*----- Check inputs to see if they are reasonable-ish -----*/
1082
1083 if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
1084
1085 if( user_func == NULL ) return NULL ;
1086
1087 if( ignore < 0 ) ignore = 0 ;
1088
1089 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
1090
1091 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum */
1092 nuse = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
1093 if( nuse < 2 ) return NULL ;
1094
1095 DSET_load( old_dset ) ; /* must be in memory before we get pointers to it */
1096
1097 kk = THD_count_databricks( old_dset->dblk ) ; /* check if it was */
1098 if( kk < DSET_NVALS(old_dset) ){ /* loaded correctly */
1099 DSET_unload( old_dset ) ;
1100 return NULL ;
1101 }
1102
1103 switch( old_datum ){ /* pointer type depends on input datum type */
1104
1105 default: /** don't know what to do **/
1106 DSET_unload( old_dset ) ;
1107 return NULL ;
1108
1109 /** create array of pointers into old dataset sub-bricks **/
1110
1111 /*--------- input is bytes ----------*/
1112 /* voxel #i at time #k is bptr[k][i] */
1113 /* for i=0..nvox-1 and k=0..nuse-1. */
1114
1115 case MRI_byte:
1116 bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
1117 if( bptr == NULL ) return NULL ;
1118 for( kk=0 ; kk < nuse ; kk++ )
1119 bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
1120 break ;
1121
1122 /*--------- input is shorts ---------*/
1123 /* voxel #i at time #k is sptr[k][i] */
1124 /* for i=0..nvox-1 and k=0..nuse-1. */
1125
1126 case MRI_short:
1127 sptr = (short **) malloc( sizeof(short *) * nuse ) ;
1128 if( sptr == NULL ) return NULL ;
1129 for( kk=0 ; kk < nuse ; kk++ )
1130 sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
1131 break ;
1132
1133 /*--------- input is floats ---------*/
1134 /* voxel #i at time #k is fptr[k][i] */
1135 /* for i=0..nvox-1 and k=0..nuse-1. */
1136
1137 case MRI_float:
1138 fptr = (float **) malloc( sizeof(float *) * nuse ) ;
1139 if( fptr == NULL ) return NULL ;
1140 for( kk=0 ; kk < nuse ; kk++ )
1141 fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
1142 break ;
1143
1144 /*--------- input is complex ---------*/
1145 /* voxel #i at time #k is cptr[k][i] */
1146 /* for i=0..nvox-1 and k=0..nuse-1. */
1147
1148 case MRI_complex:
1149 cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
1150 if( cptr == NULL ) return NULL ;
1151 for( kk=0 ; kk < nuse ; kk++ )
1152 cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
1153 break ;
1154
1155 } /* end of switch on input type */
1156
1157 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
1158
1159
1160 /*---- allocate space for 1 voxel timeseries ----*/
1161
1162 fxar = (float *) malloc( sizeof(float) * nuse ) ; /* voxel timeseries */
1163 if( fxar == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
1164
1165 /*--- get scaling factors for sub-bricks ---*/
1166
1167 fac = (float *) malloc( sizeof(float) * nuse ) ; /* factors */
1168 if( fac == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
1169
1170 use_fac = 0 ;
1171 for( kk=0 ; kk < nuse ; kk++ ){
1172 fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
1173 if( fac[kk] != 0.0 ) use_fac++ ;
1174 else fac[kk] = 1.0 ;
1175 }
1176 if( !use_fac ) ZFREEUP(fac) ;
1177
1178 /*--- setup for detrending ---*/
1179
1180 dtr = (float *) malloc( sizeof(float) * nuse ) ;
1181 if( dtr == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
1182
1183 d0fac = 1.0 / nuse ;
1184 d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
1185 for( kk=0 ; kk < nuse ; kk++ )
1186 dtr[kk] = kk - 0.5 * (nuse-1) ; /* linear trend, orthogonal to 1 */
1187
1188
1189 /*----- set up to find time at each voxel -----*/
1190
1191 tdelta = old_dset->taxis->ttdel ;
1192 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
1193 if( tdelta == 0.0 ) tdelta = 1.0 ;
1194
1195 izold = -666 ;
1196 nxy = old_dset->daxes->nxx * old_dset->daxes->nyy ;
1197
1198 /*----------------------------------------------------*/
1199 /*----- Setup has ended. Now do some real work. -----*/
1200
1201 /* start notification */
1202 #if 0
1203 user_func( 0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data ) ;
1204 #else
1205 { void (*uf)(double,double,int,float *,double,double,void *) =
1206 (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
1207 uf( 0.0l,0.0l , nvox , NULL , 0.0l,0.0l , user_data ) ;
1208 }
1209 #endif
1210
1211 /***** loop over voxels *****/
1212 for( ii=0 ; ii < nvox ; ii++ ){ /* 1 time series at a time */
1213
1214 /*** load data from input dataset, depending on type ***/
1215
1216 switch( old_datum ){
1217
1218 /*** input = bytes ***/
1219
1220 case MRI_byte:
1221 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
1222 break ;
1223
1224 /*** input = shorts ***/
1225
1226 case MRI_short:
1227 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
1228 break ;
1229
1230 /*** input = floats ***/
1231
1232 case MRI_float:
1233 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
1234 break ;
1235
1236 /*** input = complex (note we use absolute value) ***/
1237
1238 case MRI_complex:
1239 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
1240 break ;
1241
1242 } /* end of switch over input type */
1243
1244 /*** scale? ***/
1245 if( use_fac )
1246 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
1247
1248 /** compute mean and slope **/
1249
1250 x0 = x1 = 0.0 ;
1251 for( kk=0 ; kk < nuse ; kk++ ){
1252 x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
1253 }
1254
1255 x0 *= d0fac ; x1 *= d1fac ; /* factors to remove mean and trend */
1256
1257 ts_mean = x0 ;
1258 ts_slope = x1 / tdelta ;
1259
1260 /** detrend? **/
1261
1262 if( detrend )
1263 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
1264
1265 /** compute start time of this timeseries **/
1266 /* The info computed here is not being used in this version*/
1267 iz = ii / nxy ; /* which slice am I in? */
1268
1269 if( iz != izold ){ /* in a new slice? */
1270 tzero = THD_timeof( ignore ,
1271 old_dset->daxes->zzorg
1272 + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
1273 izold = iz ;
1274
1275 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
1276 }
1277
1278 /*** Send data to user function ***/
1279 #if 0
1280 user_func( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data) ;
1281 #else
1282 { void (*uf)(double,double,int,float *,double,double,void *) =
1283 (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
1284 uf( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data) ;
1285 }
1286 #endif
1287
1288
1289
1290 } /* end of outer loop over 1 voxels at a time */
1291
1292 DSET_unload( old_dset ) ;
1293
1294 /* end notification */
1295 #if 0
1296 user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data ) ;
1297 #else
1298 { void (*uf)(double,double,int,float *,double,double,void *) =
1299 (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
1300 uf( 0.0l,0.0l, 0 , NULL,0.0l,0.0l, user_data ) ;
1301 }
1302 #endif
1303
1304
1305 /*-------------- Cleanup and go home ----------------*/
1306
1307 ZFREE_WORKSPACE ;
1308 retval = 0;
1309 return &retval; /* this value is not used for now .... */
1310
1311 }
1312
1313