1 #include "SUMA_suma.h"
2 
3 static char help_msg[]= {
4 "Once you load a surface and its surface volume,,\n"
5 "its node coordinates are transformed based on the\n"
6 "surface format type and the transforms stored in\n"
7 "the surface volume. At this stage, the node coordinates\n"
8 "are in what we call RAImm DICOM where x coordinate is\n"
9 "from right (negative) to left (positive) and y coordinate\n"
10 "from anterior to posterior and z from inferior to superior\n"
11 "This RAI coordinate corresponds to the mm coordinates\n"
12 "displayed by AFNI in the top left corner of the controller\n"
13 "when you have RAI=DICOM order set (right click on coordinate\n"
14 "text are to see option. When you open the surface with the\n"
15 "same sv in SUMA and view the sv volume in AFNI, the coordinate\n"
16 "of a node on an anatomically correct surface should be close\n"
17 "to the coordinate displayed in AFNI.\n"
18 "In the output, RAImm is the coordinate just described for a \n"
19 "particular node.\n"
20 "The next coordinate in the output is called 3dfind, which stands\n"
21 "for three dimensional float index. 3dfind is a transformation \n"
22 "of the RAImm coordinates to a coordinate in the units of the\n"
23 "voxel grid. The voxel with the closest center to a location\n"
24 "at RAImm would then be at round(3dfind). In other terms, \n"
25 "RAImm is the coordinate closest to voxel  \n"
26 " V(round(3dfind[0]), round(3dfind[1]), round(3dfind[2])\n"
27 "To see index coordinates, rather than mm coordinates in \n"
28 "AFNI, set: Define Datamode --> Misc --> Voxel Coords?\n"
29 "Note that the index coordinates would be different for the\n"
30 "underlay and overlay because they are usually at different\n"
31 "resolution and/or orientation. To see the overlay coordinates\n"
32 "make sure you have 'See Overlay' turned on.\n"
33 "The last value in the output is the value from the chosen\n"
34 "sub-brick\n"
35 };
36 
usage_Surf2VolCoord_demo(SUMA_GENERIC_ARGV_PARSE * ps,int detail)37 void usage_Surf2VolCoord_demo (SUMA_GENERIC_ARGV_PARSE *ps, int detail)
38 {
39       static char FuncName[]={"usage_Surf2VolCoord_demo"};
40       char * s = NULL, *sio=NULL, *st = NULL, *sts = NULL;
41       int i;
42       s = SUMA_help_basics();
43       sio  = SUMA_help_IO_Args(ps);
44       printf ( "\n"
45 "Usage: Surf2VolCoord <-i_TYPE SURFACE> \n"
46 "                      <-grid_parent GRID_VOL> \n"
47 "                      [-grid_subbrick GSB]\n"
48 "                      [-sv SURF_VOL] \n"
49 "                      [-one_node NODE]\n"
50 "                      [-closest_nodes XYZ.1D]\n"
51 " \n"
52 "  Relates node indices to coordinates:\n"
53 "  ------------------------------------\n"
54 "  Given x y z coordinates, return the nodes closest to them.\n"
55 "  For example:\n"
56 "      Surf2VolCoord    -i SUMA/std60.lh.pial.asc \\\n"
57 "                       -i SUMA/std60.rh.pial.asc \\\n"
58 "                       -sv anat+tlrc. -qual LR   \\\n"
59 "                       -closest_nodes XYZ.1D\n"
60 "      If you are not sure you have the proper -sv, verify with SUMA:\n"
61 "        suma -i SUMA/std60.lh.pial.asc  \\\n"
62 "             -i SUMA/std60.rh.pial.asc  \\\n"
63 "             -sv anat+tlrc. -niml &\n"
64 "        afni -niml &\n"
65 "      Then press 't' in SUMA to send surfaces to AFNI.\n"
66 "\n"
67    );
68    if (detail) {
69       printf(
70 "  Mandatory Parameters:\n"
71 "     -closest_nodes XYZ.1D: A coordinate file specifying coordinates\n"
72 "                           for which the closest nodes will be found.\n"
73 /*"                           If you want to specify a coordinate on the\n"
74 "                           command line you could do:\n"
75 "                        -closest_nodes \"1D:17.2 12.4 -5.1\\'\"\n"*/
76 "                  Note: The coordinates in XYZ.1D are in RAI by default.\n"
77 "                        You can use -LPI if you need to.\n"
78 "     -closest_node 'X Y Z': An easier way to specify a single node's coords.\n"
79 "  Optional Parameters:\n"
80 "     -qual STRING: A string of characters that are used to indentify\n"
81 "                   the surface in which the closest node was found.\n"
82 "                   This is useful when you have two surfaces specified\n"
83 "                   like the left and right hemispheres for example.\n"
84 "                   In that case you can set STRING to LR if the first\n"
85 "                   surface is the left hemisphere and the second is the\n"
86 "                   right hemisphere. If you had node 12342 on the left hemi\n"
87 "                   and 7745 on the right, the output would look like this:\n"
88 "                        1342L \n"
89 "                        7745R \n"
90 "                   If qual is not set, no qualifying characters are added if \n"
91 "                   up only have one surface on the command line. \n"
92 "                   The sequence ABC... is used otherwise.\n"
93 "     -LPI: The coordinate axis direction for values in XYZ.1D are in LPI.\n"
94 "           As a result, the program will negate the sign of the X, and Y\n"
95 "           coordinates in XYZ.1D \n"
96 "     -RAI: The coordinate axis direction for values in XYZ.1D are in RAI\n"
97 "           which is the default. No transformation is applied to values in\n"
98 "           XYZ.1D\n"
99 "     -verb LEVEL: Verbosity level, default is 0\n"
100 "     -prefix PREFIX: Output results to file PREFIX (will overwrite).\n"
101 "                     Default is stdout\n"
102       );
103    }
104       printf(
105 "\n"
106 "  In Demo mode:\n"
107 "  -------------\n"
108 "  Illustrates how surface coordinates relate to voxel grid.\n"
109 "  The program outputs surface and equivalent volume coordinates\n"
110 "  for all nodes in the surface after it is aligned via its sv.\n"
111 "  The code is intended as a source code demo.\n"
112 "\n");
113    if (detail) {
114       printf(
115 "  Mandatory Parameters:\n"
116 "     -i_TYPE SURFACE: Specify input surface.\n"
117 "             You can also use -t* and -spec and -surf\n"
118 "             methods to input surfaces. See below\n"
119 "             for more details.\n"
120 "     -prefix PREFIX: Prefix of output dataset.\n"
121 "     -grid_parent GRID_VOL: Specifies the grid for the\n"
122 "                  output volume.\n"
123 "  Optional Parameters:\n"
124 "     -grid_subbrick GSB: Sub-brick from which data are taken.\n"
125 "     -one_node NODE: Output results for node NODE only.\n"
126 "\n"
127 "The output is lots of text so you're better off redirecting to a file.\n"
128       );
129    }
130    if (detail > 1) {
131        printf(
132 "%s"
133 "\n"
134 "%s"
135 "%s"
136                "\n", help_msg, sio,  s);
137       SUMA_free(s); s = NULL;
138       SUMA_free(st); st = NULL;
139       SUMA_free(sio); sio = NULL;
140    }
141    if (detail) {
142       printf("       Ziad S. Saad SSCC/NIMH/NIH saadz@mail.nih.gov     \n");
143    }
144    exit(0);
145 }
146 
SUMA_Surf2VolCoord_demo_ParseInput(char * argv[],int argc,SUMA_GENERIC_ARGV_PARSE * ps)147 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_Surf2VolCoord_demo_ParseInput(
148                      char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
149 {
150    static char FuncName[]={"SUMA_Surf2VolCoord_demo_ParseInput"};
151    SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
152    int kar, i;
153    SUMA_Boolean brk;
154    SUMA_Boolean LocalHead = NOPE;
155 
156    SUMA_ENTRY;
157 
158    Opt = SUMA_Alloc_Generic_Prog_Options_Struct();
159    Opt->obj_type = 0; /* sub-brick index */
160    Opt->NodeDbg = -1;
161    Opt->debug=0;
162    Opt->b1=0;
163    Opt->n_fvec=0;
164    kar = 1;
165    brk = NOPE;
166 	while (kar < argc) { /* loop accross command ine options */
167 		/*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
168 		if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
169 			 usage_Surf2VolCoord_demo(ps, strlen(argv[kar]) > 3 ? 2:1);
170           exit (0);
171 		}
172 
173 		SUMA_SKIP_COMMON_OPTIONS(brk, kar);
174 
175       if (!brk && (strcmp(argv[kar], "-prefix") == 0))
176       {
177          if (kar+1 >= argc)
178          {
179             fprintf (SUMA_STDERR, "need a number after -prefix \n");
180             exit (1);
181          }
182          Opt->out_vol_prefix =
183             SUMA_AfniPrefix(argv[++kar], Opt->out_vol_view, NULL,
184                             &(Opt->out_vol_exists));
185 
186          brk = YUP;
187       }
188 
189       if (!brk && (strcmp(argv[kar], "-closest_nodes") == 0))
190       {
191          if (kar+1 >= argc)
192          {
193             fprintf (SUMA_STDERR, "need a 1D file after -closest_nodes \n");
194             exit (1);
195          }
196          if (!(Opt->fvec = SUMA_Load1D_eng(argv[++kar],
197                                  &(Opt->fvec_dim), &(Opt->n_fvec), 1, 1))) {
198             SUMA_S_Errv("Failed to read %s\n", argv[kar]);
199             exit(1);
200          }
201          brk = YUP;
202       }
203 
204       if (!brk && (strcmp(argv[kar], "-closest_node") == 0))
205       {
206          if (kar+1 >= argc)
207          {
208             fprintf (SUMA_STDERR,
209                      "need a triplet in quotes after -closest_node \n");
210             exit (1);
211          }
212          Opt->fvec = (float *)SUMA_calloc(3, sizeof(float));
213          i = SUMA_StringToNum(argv[++kar], (void *)(Opt->fvec), 3, 1);
214          if (i != 3) {
215             SUMA_S_Errv("Could not get 3 values from %s\n", argv[kar]);
216             exit (1);
217          }
218          Opt->n_fvec = 1; Opt->fvec_dim = 3;
219          brk = YUP;
220       }
221 
222 
223       if (!brk && (strcmp(argv[kar], "-qual") == 0)) {
224          if (kar+1 >= argc)
225          {
226             fprintf (SUMA_STDERR, "need a string after -qual \n");
227             exit (1);
228          }
229          Opt->s = SUMA_copy_string(argv[++kar]);
230          brk = YUP;
231       }
232 
233       if (!brk && (strcmp(argv[kar], "-RAI") == 0)) {
234          Opt->b1 = 0;
235          brk = YUP;
236       }
237 
238       if (!brk && (strcmp(argv[kar], "-LPI") == 0)) {
239          Opt->b1 = 1;
240          brk = YUP;
241       }
242 
243       if (!brk && (strcmp(argv[kar], "-verb") == 0)) {
244          if (kar+1 >= argc)
245          {
246             fprintf (SUMA_STDERR, "need an integer after -verb \n");
247             exit (1);
248          }
249          Opt->debug = atoi(argv[++kar]);
250          brk = YUP;
251       }
252 
253       if (!brk && (strcmp(argv[kar], "-grid_parent") == 0))
254       {
255          if (kar+1 >= argc)
256          {
257             fprintf (SUMA_STDERR, "need a dataset after -grid_parent \n");
258             exit (1);
259          }
260          Opt->out_grid_prefix =
261             SUMA_AfniPrefix(argv[++kar], Opt->out_grid_view, NULL,
262                             &(Opt->out_grid_exists));
263          if (!SUMA_AfniExistsView(Opt->out_grid_exists, Opt->out_grid_view)) {
264             fprintf(SUMA_STDERR,
265                     "Error Surf2VolCoord_demo:\n"
266                     "Grid parent %s%s does not exist.\n",
267                     Opt->out_grid_prefix, Opt->out_grid_view);
268             exit(1);
269          }
270          brk = YUP;
271       }
272 
273       if (!brk && (strcmp(argv[kar], "-grid_subbrick") == 0))
274       {
275          if (kar+1 >= argc)
276          {
277             fprintf (SUMA_STDERR, "need a number after -grid_subbrick \n");
278             exit (1);
279          }
280          Opt->obj_type = atoi(argv[++kar]);
281          brk = YUP;
282       }
283 
284       if (!brk && (strcmp(argv[kar], "-debug") == 0))
285       {
286          if (kar+1 >= argc)
287          {
288             fprintf (SUMA_STDERR, "need a number after -debug \n");
289             exit (1);
290          }
291 
292          Opt->debug = atoi(argv[++kar]);
293          brk = YUP;
294       }
295       if (!brk && (strcmp(argv[kar], "-one_node") == 0))
296       {
297          if (kar+1 >= argc)
298          {
299             fprintf (SUMA_STDERR, "need a number after -one_node \n");
300             exit (1);
301          }
302 
303          Opt->NodeDbg = atoi(argv[++kar]);
304          brk = YUP;
305       }
306 
307       if (!brk && !ps->arg_checked[kar]) {
308 			fprintf (SUMA_STDERR,
309                   "Error Surf2VolCoord_demo:\n"
310                   "Option %s not understood. Try -help for usage\n", argv[kar]);
311 			suggest_best_prog_option(argv[0], argv[kar]);
312          exit (1);
313 		} else {
314 			brk = NOPE;
315 			kar ++;
316 		}
317    }
318 
319    SUMA_RETURN(Opt);
320 }
321 
main(int argc,char * argv[])322 int main (int argc,char *argv[])
323 {/* Main */
324    static char FuncName[]={"Surf2VolCoord_demo"};
325    SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
326    SUMA_GENERIC_ARGV_PARSE *ps=NULL;
327    SUMA_SurfSpecFile *Spec = NULL;
328    int N_Spec=0;
329    short *isin = NULL;
330    int N_in = 0, i, *minnode=NULL, *bestflag=NULL, isrf=0;
331    SUMA_SurfaceObject *SO = NULL;
332    SUMA_VOLPAR *vp = NULL;
333    THD_3dim_dataset *dset = NULL;
334    char *vpname=NULL;
335    SUMA_Boolean LocalHead = NOPE;
336    double *dvec = NULL;
337    float *tmpXYZ=NULL, *mindist=NULL;
338    int di, dj, dk, dijk, nx, ny, nxy, i0, i1, mode;
339    float fi, fj, fk;
340 
341    SUMA_STANDALONE_INIT;
342 	SUMA_mainENTRY;
343 
344    /* Allocate space for DO structure */
345 	SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
346    ps = SUMA_Parse_IO_Args(argc, argv, "-i;-t;-spec;-sv;");
347 
348    Opt = SUMA_Surf2VolCoord_demo_ParseInput (argv, argc, ps);
349    if (argc < 2) {
350       SUMA_S_Err("Too few options");
351       usage_Surf2VolCoord_demo(ps, 0);
352       exit (1);
353    }
354 
355 
356    if (Opt->debug > 2) LocalHead = YUP;
357    if (Opt->n_fvec) mode = 1;
358    else mode = 0;
359 
360    /* some checks ...*/
361    if (!Opt->out_vol_prefix) {
362       if (mode==0) {
363          Opt->out_vol_prefix =
364             SUMA_AfniPrefix("Surf2VolCoord_demo",
365                             NULL, NULL, &(Opt->out_vol_exists));
366          strncpy(Opt->out_vol_view, Opt->out_grid_view, SUMA_VIEW_LENGTH);
367       }
368    } else {
369       if (SUMA_AfniExistsView(Opt->out_vol_exists, Opt->out_vol_view)) {
370          fprintf(SUMA_STDERR,
371                   "Error Surf2VolCoord_demo:\nOutput volume %s%s exists.\n",
372                   Opt->out_vol_prefix, Opt->out_vol_view);
373          exit(1);
374       }
375    }
376 
377    if (Opt->out_grid_prefix) {
378       if (!SUMA_AfniExistsView(Opt->out_grid_exists, Opt->out_grid_view)) {
379          fprintf(SUMA_STDERR, "Error Surf2VolCoord_demo:\n"
380                               "Grid parent %s%s does not exist.\n",
381                               Opt->out_grid_prefix, Opt->out_grid_view);
382          exit(1);
383       }
384    }
385 
386    /* check on inputs */
387    if ((ps->s_N_surfnames && ps->i_N_surfnames) ||
388        (ps->s_N_surfnames &&  ps->t_N_surfnames) ||
389        (ps->i_N_surfnames &&  ps->t_N_surfnames)) {
390       SUMA_S_Err("Multiple surface input methods used."
391                  "Only one method allowed.");
392       exit(1);
393    }
394 
395 
396    Spec = SUMA_IO_args_2_spec(ps, &N_Spec);
397    if (N_Spec == 0) {
398       SUMA_S_Err("No surfaces found.");
399       exit(1);
400    }
401 
402    if (N_Spec != 1) {
403       SUMA_S_Err("Multiple spec at input.");
404       exit(1);
405    }
406 
407    if (mode == 0 && Spec->N_Surfs!=1) {
408       SUMA_S_Err("Only one surface allowed in demo mode");
409       exit(1);
410    }
411    if (Opt->s && strlen(Opt->s) < Spec->N_Surfs) {
412       SUMA_S_Errv("Have %d surfaces but qualifier string %s is %d characters.\n",
413                   Spec->N_Surfs, Opt->s, (int)strlen(Opt->s));
414       exit(1);
415    }
416 
417    if (mode == 1 && Opt->b1) {/* go from LPI to RAI */
418       for (i=0; i<Opt->n_fvec; ++i) {
419          Opt->fvec[Opt->fvec_dim*i] = -Opt->fvec[Opt->fvec_dim*i];
420          Opt->fvec[Opt->fvec_dim*i+1] = -Opt->fvec[Opt->fvec_dim*i+1];
421       }
422    }
423 
424    for (isrf=0; isrf < Spec->N_Surfs; ++isrf) {
425       SO = SUMA_Load_Spec_Surf(Spec, isrf, ps->sv[0], 0);
426       if (!SO) {
427             fprintf (SUMA_STDERR,"Error %s:\n"
428                                  "Failed to find surface\n"
429                                  "in spec file. \n",
430                                  FuncName );
431             exit(1);
432 
433       }
434       if (Opt->NodeDbg >= 0 && Opt->NodeDbg >= SO->N_Node) {
435          fprintf (SUMA_STDERR,"Error %s:\n"
436                                  "Node index %d is >= SO->N_Node (%d)\n"
437                                  "\n",
438                                  FuncName, Opt->NodeDbg, SO->N_Node );
439             exit(1);
440       }
441 
442 
443       if (mode == 1) {
444          double dd=0.0;
445          int n;
446          if (Opt->fvec_dim != SO->NodeDim) {
447             SUMA_S_Errv("Have %d dims for xyz, and %d for surf.\n",
448                   Opt->fvec_dim, SO->NodeDim);
449             exit (1);
450          }
451          if (!mindist) mindist =
452                (float *)SUMA_calloc(Opt->n_fvec,sizeof(float));
453          if (!minnode) minnode =
454                (int *)SUMA_calloc(Opt->n_fvec,sizeof(int));
455          if (!bestflag) bestflag =
456                (int *)SUMA_calloc(Opt->n_fvec,sizeof(int));
457          for (i=0; i<Opt->n_fvec; ++i) {
458             if (isrf==0) mindist[i] = 156779800000.0;
459             for (n=0; n<SO->N_Node; ++n) {
460                SUMA_SEG_NORM( (SO->NodeList+SO->NodeDim*n),
461                               (Opt->fvec+Opt->fvec_dim*i), dd);
462                if (dd < mindist[i]) {
463                   mindist[i] = dd; minnode[i] = n; bestflag[i]=isrf;
464                }
465             }
466          }
467       } else {
468          /***********          The demo mode *****************/
469          /* By now SO is the surface object whose coordinates have transformed
470          so that it is in register with the surface volume specifed on command
471          line.
472          */
473 
474          /* Now let us read the volume from which you would be accessing values.
475          That volume should be in alignment with -sv but not at the same
476          resolution.
477          I assume that this volume contains one sub-brick for simplicity. If no
478          such volume is specified then we'll use the -sv volume for that */
479 
480          if (Opt->out_grid_prefix) {
481             vpname = SUMA_append_string(Opt->out_grid_prefix,
482                                         Opt->out_grid_view);
483             vp = SUMA_VolPar_Attr(vpname);
484             dset = THD_open_dataset( vpname );
485          } else {
486             vp = SO->VolPar;
487             if (!vp) {
488                fprintf (SUMA_STDERR,"Error %s:\n"
489                                        "Need a grid parent.\n",
490                                        FuncName );
491                   exit(1);
492             }
493             vpname = SUMA_copy_string(ps->sv[0]);
494             if (!SUMA_AfniView(ps->sv[0], Opt->out_grid_view)) {
495                fprintf (SUMA_STDERR,"Error %s:\n"
496                                        "Failed to get view!!!\n",
497                                        FuncName );
498                   exit(1);
499 
500             }
501             dset = THD_open_dataset( vpname );
502          }
503          /* load .BRIK into memory */
504          DSET_load(dset);
505          if (LocalHead) {
506             fprintf(SUMA_STDERR,"%s: Using %s for grid\n", FuncName, vpname);
507          }
508          if (Opt->obj_type >= DSET_NVALS(dset)) {
509             fprintf (SUMA_STDERR,"Error %s:\n"
510                                  "Grid dset has a total of %d sub-bricks\n"
511                                  " user specified sb#%d!!!\n",
512                                  FuncName, DSET_NVALS(dset), Opt->obj_type );
513                   exit(1);
514          }
515          /* transform the surface's coordinates from RAImm to 3dfind
516             RAI coordinates are the mm coordinates you see displayed in AFNI's
517                top left corner
518             3dfind are the equivalent coordinates expressed in index coordinates
519                on the grid volume
520             Say you have a node at 13.4, 23, -54 mm (RAI), its 3dfind coordinate
521                might be: 3.2 , 56.1,   124.8 which would be closest to voxel with
522                ijk indices 3, 56, 125 in the volume defined by grid parent
523                (dset or vp).
524             To get that voxel's value, you just access element 3,56, 125.
525             Note that in this example, I am getting voxel values for all  coords
526             in vector tmpXYZ and those happen to be all the nodes of the surface.
527             But you could put whatever coordinates or subset of coordinates you
528             want in there. Just be sure to change SO->N_Node to reflect
529             the number of triplets in tmpXYZ
530          */
531 
532          /* copy surface coordinates to preserve them, we're going to ijk land */
533          tmpXYZ = (float *)SUMA_malloc(SO->N_Node * 3 * sizeof(float));
534          if (!tmpXYZ) {
535             SUMA_SL_Crit("Faile to allocate");
536             exit(1);
537          }
538          memcpy ((void*)tmpXYZ, (void *)SO->NodeList,
539                   SO->N_Node * 3 * sizeof(float));
540          if (!SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, vp)) {
541             SUMA_SL_Err("Failed to effectuate coordinate transform.");
542             SUMA_free(tmpXYZ); tmpXYZ = NULL;
543             exit(1);
544          }
545 
546          /* Now, let us loop through all the coordinates of interest and see what
547          voxels they fall in and what values are at those voxels. */
548 
549          /* first let us load the data from one sub-brick into a double vector
550             (recall that data can be stored in variety of precisions on disk). */
551          dvec = (double *)SUMA_malloc(sizeof(double) * DSET_NVOX(dset));
552          if (!dvec) {
553             SUMA_S_Errv("Failed to allocate for %d dvec.\nOh misery.\n",
554                         DSET_NVOX(dset));
555             exit(1);
556          }
557          EDIT_coerce_scale_type( DSET_NVOX(dset) ,
558                                  DSET_BRICK_FACTOR(dset,Opt->obj_type) ,
559                                  DSET_BRICK_TYPE(dset,Opt->obj_type),
560                                  DSET_ARRAY(dset, Opt->obj_type) ,   /* input  */
561                                  MRI_double               , dvec  ) ;/* output */
562 
563          nx = DSET_NX(dset); ny = DSET_NY(dset); nxy = nx * ny;
564          if (Opt->NodeDbg >= 0) { i0 = Opt->NodeDbg; i1 = Opt->NodeDbg+1; }
565          else { i0 = 0; i1 = SO->N_Node;};
566          for (i=i0; i<i1; ++i) {
567             fi = tmpXYZ[3*i  ]; di = SUMA_ROUND(fi);
568             fj = tmpXYZ[3*i+1]; dj = SUMA_ROUND(fj);
569             fk = tmpXYZ[3*i+2]; dk = SUMA_ROUND(fk);
570             dijk = SUMA_3D_2_1D_index(di, dj, dk, nx, nxy);
571             fprintf(SUMA_STDOUT,
572                      "Node Index %d: RAImm %.3f %.3f %.3f : "
573                      "3dfind %.1f %.1f %.1f : 3dind %d %d %d : Val %f\n",
574                       i,
575                    SO->NodeList[3*i], SO->NodeList[3*i+1], SO->NodeList[3*i+2],
576                       fi, fj, fk, di, dj, dk,
577                       dvec[dijk]);
578 
579          }
580          SUMA_free(dvec); dvec = NULL;
581          SUMA_free(tmpXYZ); tmpXYZ = NULL;
582          /* no need for data in input volume anymore */
583          PURGE_DSET(dset);
584       }
585 
586       if (vpname) SUMA_free(vpname); vpname = NULL;
587       if (dset) { DSET_delete(dset); dset = NULL; }
588       if (vp != SO->VolPar) SUMA_Free_VolPar(vp); vp = NULL;
589       if (SO) SUMA_Free_Surface_Object(SO); SO = NULL;
590    }
591 
592    if (Opt->n_fvec) {
593       FILE *fout=NULL;
594       if (Opt->out_vol_prefix) {
595          if (!(fout = fopen(Opt->out_vol_prefix,"w"))) {
596             SUMA_S_Errv("Failed to open %s for output",
597                         Opt->out_vol_prefix);
598             exit(1);
599          }
600       } else {
601          fout = SUMA_STDOUT;
602       }
603       if (Opt->debug) fprintf(fout,"#Closest_Node Distance_to_Coord\n");
604       for (i=0; i<Opt->n_fvec; ++i) {
605          fprintf(fout,"%d%c %f\n",
606                   minnode[i],
607                   Opt->s ? Opt->s[bestflag[i]]:'A'+bestflag[i],
608                   mindist[i]);
609       }
610       if (fout != SUMA_STDOUT) fclose(fout); fout=NULL;
611    }
612 
613 
614    if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
615    if (N_Spec) {
616       int k=0;
617       for (k=0; k<N_Spec; ++k) {
618          if (!SUMA_FreeSpecFields(&(Spec[k]))) {
619             SUMA_S_Err("Failed to free spec fields"); }
620       }
621       SUMA_free(Spec); Spec = NULL; N_Spec = 0;
622    }
623    if (mindist) SUMA_free(mindist); mindist = NULL;
624    if (minnode) SUMA_free(minnode); minnode = NULL;
625    if (bestflag) SUMA_free(bestflag); bestflag = NULL;
626    if (Opt) Opt = SUMA_Free_Generic_Prog_Options_Struct(Opt);
627    if (!SUMA_Free_CommonFields(SUMAg_CF))
628       SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
629    exit(0);
630 
631 }
632