#include "SUMA_suma.h"


/*#define  DO_SCALE_RANGE   *//*!< scale node coordinates to 0 <--> 100. DO NOT USE IT, OBSOLETE*/
#ifndef DO_SCALE_RANGE
   #define DO_SCALE 319.7   /*!< scale node coordinates by specified factor. Useful for tesscon coordinate system in iv files*/
#endif

SUMA_SurfaceViewer *SUMAg_cSV; /*!< Global pointer to current Surface Viewer structure*/
SUMA_SurfaceViewer *SUMAg_SVv; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
                                    SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
SUMA_DO *SUMAg_DOv;   /*!< Global pointer to Displayable Object structure vector*/
int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
SUMA_CommonFields *SUMAg_CF; /*!< Global pointer to structure containing info common to all viewers */

void usage_SUMA_SurfaceMetrics ()
   {
      static char FuncName[]={"usage_SUMA_SurfaceMetrics"};
      char * s = NULL;
      s = SUMA_help_basics();
      printf ( "\n"
               "Usage: SurfaceMetrics <-Metric1> [[-Metric2] ...] \n"
               "                  <-spec SpecFile> <-surf_A insurf> \n"
               "                  [<-sv SurfaceVolume [VolParam for sf surfaces]>]\n"
               "                  [-tlrc] [<-prefix prefix>]\n"
               "\n"
               "Outputs information about a surface's mesh\n"
               "\n"
               "   -Metric1: Replace -Metric1 with the following:\n"
               "      -vol: calculates the volume of a surface.\n"
               "            Volume unit is the cube of your surface's\n"
               "            coordinates unit, obviously.\n"
               "            Volume's sign depends on the orientation\n"
               "            of the surface's mesh.\n" 
               "            Make sure your surface is a closed one\n"
               "            and that winding is consistent.\n"
               "            Use SurfQual to check the surface.\n"
               "            If your surface's mesh has problems,\n"
               "            the result is incorrect. \n"
               "            Volume is calculated using Gauss's theorem,\n"
               "            see [Hughes, S.W. et al. 'Application of a new \n"
               "            discreet form of Gauss's theorem for measuring \n"
               "            volume' in Phys. Med. Biol. 1996].\n"
               "      -conv: output surface convexity at each node.\n"
               "         Output file is prefix.conv. Results in two columns:\n"
               "         Col.0: Node Index\n"
               "         Col.1: Convexity\n"
               "         This is the measure used to shade sulci and gyri in SUMA.\n"
               "         C[i] = Sum(dj/dij) over all neighbors j of i\n"
               "         dj is the distance of neighboring node j to the tangent plane at i\n"
               "         dij is the length of the segment ij\n"
               "      -closest_node XYZ_LIST.1D: Find the closest node on the surface\n"
               "                              to each XYZ triplet in XYZ_LIST.1D\n"
               "                              Note that it is assumed that the XYZ\n"
               "                              coordinates are in RAI (DICOM) per AFNI's\n"
               "                              coordinate convention. For correspondence\n"
               "                              with coordinates observed in SUMA and AFNI\n"
               "                              be sure to use the proper -sv parameter for\n"
               "                              the surface and XYZ coordinates in question.\n"
               "         Output file is prefix.closest.1D. Results in 8 columns:\n"
               "         Col.0: Index of closest node.\n"
               "         Col.1: Distance of closest node to XYZ reference point.\n"
               "         Col.2..4: XYZ of reference point (same as XYZ_LIST.1D, copied \n"
               "                   here for clarity).\n"
               "         Col.5..7: XYZ of closest node (after proper surface coordinate\n"
               "                   transformation, including SurfaceVolume transform.\n"
               "      -area: output area of each triangle. \n"
               "         Output file is prefix.area. Results in two columns:\n"
               "         Col.0: Triangle Index\n"
               "         Col.1: Triangle Area\n"
               "      -curv: output curvature at each node.\n"
               "         Output file is prefix.curv. Results in nine columns:\n"
               "         Col.0: Node Index\n"
               "         Col.1-3: vector of 1st principal direction of surface\n"
               "         Col.4-6: vector of 2nd principal direction of surface\n"
               "         Col.7: Curvature along T1\n"
               "         Col.8: Curvature along T2\n"
               "         Curvature algorithm by G. Taubin from: \n"
               "         'Estimating the tensor of curvature of surface \n"
               "         from a polyhedral approximation.'\n"
               "      -edges: outputs info on each edge. \n"
               "         Output file is prefix.edges. Results in five columns:\n"
               "         Col.0: Edge Index (into a SUMA structure).\n"
               "         Col.1: Index of the first node forming the edge\n"
               "         Col.2: Index of the second node forming the edge\n"
               "         Col.3: Number of triangles containing edge\n"
               "         Col.4: Length of edge.\n"
               "      -node_normals: Outputs segments along node normals.\n"
               "                     Segments begin at node and have a default\n"
               "                     magnitude of 1. See option 'Alt+Ctrl+s' in \n"
               "                     SUMA for visualization.\n"
               "      -face_normals: Outputs segments along triangle normals.\n"
               "                     Segments begin at centroid of triangles and \n"
               "                     have a default magnitude of 1. See option \n"
               "                     'Alt+Ctrl+s' in SUMA for visualization.\n"
               "      -normals_scale SCALE: Scale the normals by SCALE (1.0 default)\n"
               "                     For use with options -node_normals and -face_normals\n"
               "      -coords: Output coords of each node after any transformation \n"
               "         that is normally carried out by SUMA on such a surface.\n"
               "         Col. 0: Node Index\n"
               "         Col. 1: X\n"
               "         Col. 2: Y\n"
               "         Col. 3: Z\n"     
               "      -sph_coords: Output spherical coords of each node.\n"
               "      -sph_coords_center x y z: Shift each node by  x y z\n"
               "                                before calculating spherical\n"
               "                                coordinates. Default is the\n"
               "                                center of the surface.\n"
               "          Both sph_coords options output the following:\n"
               "          Col. 0: Node Index\n"
               "          Col. 1: R (radius)\n"
               "          Col. 2: T (azimuth)\n"
               "          Col. 3: P (elevation)\n"
               "      -boundary_nodes: Output nodes that form a boundary of a surface\n"
               "                   i.e. they form edges that belong to one and only\n"
               "                   one triangle.\n"
               "      -internal_nodes: Output nodes that are not a boundary.\n"
               "                   i.e. they form edges that belong to more than\n"
               "                   one triangle.\n"
               "\n"
               "      You can use any or all of these metrics simultaneously.\n"
               "\n"
               "   -spec SpecFile: Name of specfile containing surface of interest.\n"
               "                   If the surface does not have a spec file, use the \n"
               "                   program quickspec to create one.\n"
               "   -surf_A insurf: Name of surface of interest. \n"
               "                   NOTE: i_TYPE inSurf option is not supported for this program.\n"
               "\n"
               "   -sv SurfaceVolume [VolParam for sf surfaces]: Specify a surface volume\n"
               "                   for surface alignment. See ConvertSurface -help for more info.\n"
               "\n"
               "   -tlrc: Apply Talairach transform to surface.\n"
               "                   See ConvertSurface -help for more info.\n"
               "\n"
               "   -prefix prefix: Use prefix for output files. (default is prefix of inSurf)\n"
               "%s"
               "\n", s);
      SUMA_free(s); s = NULL;
      s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
      printf ( "       Ziad S. Saad SSCC/NIMH/NIH saadz@mail.nih.gov \n"
               "       Mon May 19 15:41:12 EDT 2003\n"
               "\n");   
   }
#define SURFACEMETRICS_MAX_SURF 10

int main (int argc,char *argv[])
{/* Main */
   static char FuncName[]={"Main_SUMA_SurfaceMetrics"};
   char  *OutName=NULL, *OutPrefix = NULL, *if_name = NULL, 
         *if_name2 = NULL, *sv_name = NULL, *vp_name = NULL,
         *tlrc_name = NULL;
   float *Cx = NULL, sph_center[3], NormScale;
   SUMA_STRING *MetricList = NULL;
   int i, n1, n2, n1_3, n2_3, kar, nt, SO_read;
   double edgeL2;
   FILE *fout=NULL;
   SUMA_SO_File_Type iType = SUMA_FT_NOT_SPECIFIED;
   SUMA_SurfaceObject *SO = NULL;   
   SUMA_SFname *SF_name = NULL;
   void *SO_name = NULL;   
   SUMA_SurfSpecFile Spec;
   THD_warp *warp=NULL ;
   THD_3dim_dataset *aset=NULL;
   char *surf_names[SURFACEMETRICS_MAX_SURF];
   char *spec_file, *histnote;
   char *closest_to_xyz = NULL;
   int insurf_method = 0, N_surf = 0, ind = 0;
   SUMA_Boolean   brk, Do_tlrc, Do_conv, Do_curv, 
                  Do_area, Do_edges, Do_vol, Do_sph, NewCent, Do_cord, Do_TriNorm, 
                  Do_NodeNorm, Do_en, Do_in, LocalHead = NOPE;  
   
	SUMA_STANDALONE_INIT;
   SUMA_mainENTRY;
   
   
	/* Allocate space for DO structure */
	SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);

   if (argc < 4)
       {
          SUMA_S_Err("Too few parameters");
          usage_SUMA_SurfaceMetrics ();
          exit (1);
       }
   
   MetricList = SUMA_StringAppend (NULL, NULL);
   kar = 1;
	brk = NOPE;
   Do_cord = NOPE;
   Do_sph = NOPE;
   Do_vol = NOPE;
   Do_tlrc = NOPE;
   Do_conv = NOPE;
   Do_area = NOPE;
   Do_curv = NOPE;
   Do_edges = NOPE;
   Do_TriNorm = NOPE;
   Do_NodeNorm = NOPE;
   Do_en = NOPE;
   Do_in = NOPE;
   closest_to_xyz = NULL;
   NormScale = 5.0;
   NewCent = NOPE;
   OutPrefix = NULL;
   for (i=0; i<SURFACEMETRICS_MAX_SURF; ++i) { surf_names[i] = NULL; }   
   spec_file = NULL;

	while (kar < argc) { /* loop accross command ine options */
		/* fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName); */
		if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
			 usage_SUMA_SurfaceMetrics();
          exit (1);
		}
		
      SUMA_SKIP_COMMON_OPTIONS(brk, kar);
      
		if (!brk && (strcmp(argv[kar], "-i_fs") == 0)) {
         SUMA_SL_Err("Option -i_fs is not supported for this program.\nUse -spec and -surf_A instead.\n");
         exit(1);
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -i_fs ");
				exit (1);
			}
			if_name = argv[kar];
         iType = SUMA_FREE_SURFER;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-i_sf") == 0)) {
         SUMA_SL_Err("Option -i_sf is not supported for this program.\nUse -spec and -surf_A instead.\n");
         exit(1);
         kar ++;
			if (kar+1 >= argc)  {
		  		fprintf (SUMA_STDERR, "need 2 argument after -i_sf");
				exit (1);
			}
			if_name = argv[kar]; kar ++;
         if_name2 = argv[kar];
         iType = SUMA_SUREFIT;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-i_vec") == 0)) {
         SUMA_SL_Err("Option -i_vec is not supported for this program.\nUse -spec and -surf_A instead.\n");
         exit(1);
         kar ++;
			if (kar+1 >= argc)  {
		  		fprintf (SUMA_STDERR, "need 2 argument after -i_vec");
				exit (1);
			}
			if_name = argv[kar]; kar ++;
         if_name2 = argv[kar];
         iType = SUMA_VEC;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-i_ply") == 0)) {
         SUMA_SL_Err("Option -i_ply is not supported for this program.\nUse -spec and -surf_A instead.\n");
         exit(1);
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -i_ply ");
				exit (1);
			}
			if_name = argv[kar];
         iType = SUMA_PLY;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-spec") == 0)) {
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -spec \n");
				exit (1);
			}
			spec_file = argv[kar];
         if (!insurf_method) insurf_method = 2;
         else {
            fprintf (SUMA_STDERR, "already specified spec file.\n");
            exit(1);
         }
			brk = YUP;
		}
      
      if (!brk && (strncmp(argv[kar], "-surf_", 6) == 0)) {
			if (kar + 1>= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -surf_X SURF_NAME \n");
				exit (1);
			}
			ind = argv[kar][6] - 'A';
         if (ind < 0 || ind >= SURFACEMETRICS_MAX_SURF) {
            fprintf (SUMA_STDERR, "-surf_X SURF_NAME option is out of range,\n"
                                  "   only surf_A allowed.\n");
				exit (1);
         }
         kar ++;
         surf_names[ind] = argv[kar];
         N_surf = ind+1;
			if (insurf_method != 2) {
            fprintf (SUMA_STDERR, "-surf_X SURF_NAME option must be used with -spec option.\n");
            exit(1);
         }
         brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-sph_coords_center") == 0)) {
         kar ++;
			if (kar+2 >= argc)  {
		  		fprintf (SUMA_STDERR, "need 3 arguments after -sph_coords_center \n");
				exit (1);
			}
			sph_center[0] = atof(argv[kar]); kar ++;
         sph_center[1] = atof(argv[kar]); kar ++;
         sph_center[2] = atof(argv[kar]);
         NewCent = YUP;
         Do_sph = YUP;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-sph_coords") == 0)) {
         Do_sph = YUP;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-coords") == 0)) {
         Do_cord = YUP;
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-prefix") == 0)) {
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -prefix ");
				exit (1);
			}
			OutPrefix = SUMA_copy_string(argv[kar]);
			brk = YUP;
		}
      
      if (!brk && (strcmp(argv[kar], "-sv") == 0)) {
         if (0 && iType == SUMA_FT_NOT_SPECIFIED) {/* iType input no longer allowed */
            fprintf (SUMA_STDERR, " -sv option must be preceeded by -i_TYPE option.");
            exit(1);
         }
         kar ++;
			if (iType == SUMA_SUREFIT) {
            if (kar+1 >= argc)  {
		  		   fprintf (SUMA_STDERR, "need 2 argument after -sv (SurfaceVolume and VolumeParent)");
				   exit (1);
			   }
            sv_name = argv[kar]; kar ++;
            vp_name = argv[kar];
         } else {
            if (kar >= argc)  {
		  		   fprintf (SUMA_STDERR, "need argument after -sv ");
				   exit (1);
			   }
			   sv_name = argv[kar];
         }
			brk = YUP;
		}

      if (!brk && (strcmp(argv[kar], "-tlrc") == 0)) {
         Do_tlrc = YUP;
         brk = YUP;
      }
      if (!brk && (strcmp(argv[kar], "-node_normals") == 0)) {
         Do_NodeNorm = YUP;
         brk = YUP;
      }
      if (!brk && (strcmp(argv[kar], "-face_normals") == 0)) {
         Do_TriNorm = YUP;
         brk = YUP;
      }
      if (!brk && (strcmp(argv[kar], "-normals_scale") == 0)) {
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -normals_scale ");
				exit (1);
			}
         NormScale = atof(argv[kar]);
         brk = YUP;
      }
      if (!brk && (strcmp(argv[kar], "-closest_node") == 0)) {
         kar ++;
			if (kar >= argc)  {
		  		fprintf (SUMA_STDERR, "need argument after -closest_node ");
				exit (1);
			}
         closest_to_xyz = argv[kar];
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-conv") == 0)) {
         Do_conv = YUP;
         MetricList = SUMA_StringAppend (MetricList, "Convexity "); 
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-area") == 0)) {
         Do_area = YUP;
         MetricList = SUMA_StringAppend (MetricList, "PolyArea "); 
         brk = YUP;
      }

      if (!brk && (strcmp(argv[kar], "-curv") == 0)) {
         Do_curv = YUP;
         MetricList = SUMA_StringAppend (MetricList, "Curvature "); 
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-edges") == 0)) {
         Do_edges = YUP;
         MetricList = SUMA_StringAppend (MetricList, "EdgeList "); 
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-vol") == 0)) {
         Do_vol = YUP;
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-boundary_nodes") == 0)) {
         Do_en = YUP;
         brk = YUP;
      }
      
      if (!brk && (strcmp(argv[kar], "-internal_nodes") == 0)) {
         Do_in = YUP;
         brk = YUP;
      }
      
      if (!brk) {
			fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
			exit (1);
		} else {	
			brk = NOPE;
			kar ++;
		}
   }
   
   if (N_surf < 1) {
      SUMA_SL_Err("No surface specified.");
      exit(1);
   }

   /* clean MetricList */
   MetricList = SUMA_StringAppend (MetricList, NULL); 
   
   /* sanity checks */
   if (!strlen(MetricList->s) && !Do_vol && !Do_sph && !Do_cord && !Do_TriNorm && !Do_NodeNorm && !Do_en && !Do_in && !closest_to_xyz) {
      SUMA_S_Err("No Metrics specified.\nNothing to do.\n");
      exit(1);
   }
   
   if (0 && sv_name) { /* stupid check for volumes... */
      if (!SUMA_filexists(sv_name)) {
         fprintf (SUMA_STDERR,"Error %s: %s not found.\n", FuncName, sv_name);
         exit(1);
      }
   }
   
   if (Do_tlrc && !sv_name) {
      fprintf (SUMA_STDERR,"Error %s: -tlrc must be used with -sv option.\n", FuncName);
      exit(1);
   }
   
   if (vp_name) {
      if (!SUMA_filexists(vp_name)) {
         fprintf (SUMA_STDERR,"Error %s: %s not found.\n", FuncName, vp_name);
         exit(1);
      }
   }

   
   /* read all surfaces */
   if (!SUMA_AllocSpecFields(&Spec)) { SUMA_S_Err("Error initing"); exit(1); }
   if (!SUMA_Read_SpecFile (spec_file, &Spec)) {
		fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
		exit(1);
	}
   SO_read = SUMA_spec_select_surfs(&Spec, surf_names, SURFACEMETRICS_MAX_SURF, 0);
   if ( SO_read != N_surf )
   {
	   if (SO_read >=0 )
         fprintf(SUMA_STDERR,"Error %s:\nFound %d surfaces, expected %d.\n", FuncName,  SO_read, N_surf);
      exit(1);
   }
   
   /* now read into SUMAg_DOv */
   if (!SUMA_LoadSpec_eng(&Spec, SUMAg_DOv, &SUMAg_N_DOv, sv_name, 0, SUMAg_CF->DsetList) ) {
	   fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_LoadSpec_eng\n", FuncName);
      exit(1);
   }   SO = SUMA_find_named_SOp_inDOv(surf_names[0], SUMAg_DOv, SUMAg_N_DOv);
   if (!SO) {
      fprintf (SUMA_STDERR,"Error %s: Failed to read input surface.\n", FuncName);
      exit (1);
   }
   
   if (Do_tlrc) {
      fprintf (SUMA_STDOUT,"Performing talairach transform...\n");

      /* form the tlrc version of the surface volume */
      tlrc_name = (char *) SUMA_calloc (strlen(SO->VolPar->dirname)+strlen(SO->VolPar->prefix)+60, sizeof(char));
      sprintf (tlrc_name, "%s%s+tlrc.HEAD", SO->VolPar->dirname, SO->VolPar->prefix);
      if (!SUMA_filexists(tlrc_name)) {
         fprintf (SUMA_STDERR,"Error %s: %s not found.\n", FuncName, tlrc_name);
         exit(1);
      }
      
      /* read the tlrc header */
      aset = THD_open_dataset(tlrc_name) ;
      if( !ISVALID_DSET(aset) ){
         fprintf (SUMA_STDERR,"Error %s: %s is not a valid data set.\n", FuncName, tlrc_name) ;
         exit(1);
      }
      if( aset->warp == NULL ){
         fprintf (SUMA_STDERR,"Error %s: tlrc_name does not contain a talairach transform.\n", FuncName);
         exit(1);
      }
      
      warp = aset->warp ;
      
      /* now warp the coordinates, one node at a time */
      if (!SUMA_AFNI_forward_warp_xyz(warp, SO->NodeList, SO->N_Node)) {
         fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_AFNI_forward_warp_xyz.\n", FuncName);
         exit(1);
      }

      
      
   }
   
   /* check on normals */
   if ((Do_TriNorm || Do_NodeNorm) && (!SO->NodeNormList || !SO->FaceNormList)) {
      SUMA_RECOMPUTE_NORMALS(SO);
   }
   
   /* create the surface label*/
   SO->Label = SUMA_SurfaceFileName (SO, NOPE);
   if (!SO->Label) {
      SUMA_S_Err("Failed to create Label");
      exit(1);
   }

   if (LocalHead) SUMA_Print_Surface_Object (SO, stderr);
   
   /* Now do the deed */
   SUMA_LH (MetricList->s);
   if (strlen(MetricList->s)) {
      if (!SUMA_SurfaceMetrics (SO, MetricList->s, NULL)) {
         fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_SurfaceMetrics.\n", FuncName);
         exit(1);
      }
   }

   SUMA_LH ("Done with Metrics");
   
   /* output time */
   if (!OutPrefix) {
      OutPrefix = SUMA_copy_string(SO->Label);
   }
   
   OutName = (char*) SUMA_malloc((strlen(OutPrefix) + 30) * sizeof(char));
   histnote = SUMA_HistString (NULL, argc, argv, NULL);
   
   if (Do_sph) {
      double *sph=NULL;
      sprintf(OutName, "%s.sphcoord.1D.dset", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Edge output file exists.\nWill not overwrite.");
         exit(1);
      }
      if (NewCent) ;
      else if (SO->Center) { SUMA_COPY_VEC(SO->Center, sph_center, 3, float, float); }
      else {
         SUMA_SL_Err("SO has no center");
         exit(1);
      }
      sph = SUMA_Cart2Sph(SO->NodeList, SO->N_Node, sph_center);

      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#Spherical coords, \n");
      fprintf (fout,"#  cartesian coords shifted by [%f %f %f] prior to xform\n", sph_center[0], sph_center[1], sph_center[2]);
      fprintf (fout,"#nI = Node Index\n");
      fprintf (fout,"#r  = Rho (radius)\n");
      fprintf (fout,"#t  = theta(azimuth)\n");
      fprintf (fout,"#p  = phi(elevation)\n");
      fprintf (fout,"#nI\tr\tt\tp\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i < SO->N_Node; ++i) {
            fprintf (fout,"%d\t%f\t%f\t%f\n",
                  i, sph[3*i], sph[3*i+1],sph[3*i+2]);
      }
      
      fclose(fout); fout = NULL;
      
      if (sph) SUMA_free(sph); sph = NULL;
   }  
   
   if (Do_NodeNorm) {
      float norm[3];
      sprintf(OutName, "%s.NodeNormSeg.1D", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Node normals output file exists.\nWill not overwrite.");
         exit(1);
      }
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#Node normals.\n");
      fprintf (fout,"#  Segments from node along the direction of the normal (of magnitude %f)\n", NormScale);
      fprintf (fout,"#  1st three columns are node's coordinates\n");
      fprintf (fout,"#  2nd three columns are the coordinate of a second point along the normal\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i<SO->N_Node; ++i) {
         norm[0] = SO->NodeNormList[3*i]; norm[1] = SO->NodeNormList[3*i+1]; norm[2] = SO->NodeNormList[3*i+2];
         /* normal is expected to be normalized...*/
         norm[0] *= NormScale; norm[1] *= NormScale; norm[2] *= NormScale;
         fprintf (fout,"%f %f %f \t%f %f %f\n", 
            SO->NodeList[3*i], SO->NodeList[3*i+1], SO->NodeList[3*i+2], 
            SO->NodeList[3*i]+norm[0], SO->NodeList[3*i+1]+norm[1], SO->NodeList[3*i+2]+norm[2]);            
      }
      
      fclose(fout); fout = NULL;
      
   }
   
   if (Do_TriNorm) {
      float tc[3], norm[3];
      int n1, n2, n3;
      if (SO->FaceSetDim != 3) {
         SUMA_S_Err("Triangular meshes only please.");
         exit(1);
      }
      sprintf(OutName, "%s.TriNormSeg.1D", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Triangle normals output file exists.\nWill not overwrite.");
         exit(1);
      }
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#Triangle normals.\n");
      fprintf (fout,"#  Segments from centroid of triangle along the direction of the normal (of magnitude %f)\n", NormScale);
      fprintf (fout,"#  1st three columns are centroid's coordinates\n");
      fprintf (fout,"#  2nd three columns are the coordinate of a second point along the normal\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i<SO->N_FaceSet; ++i) {
         n1 = SO->FaceSetList[3*i]; n2 = SO->FaceSetList[3*i+1]; n3 = SO->FaceSetList[3*i+2];
         /* coordinate of centroid */
         tc[0] = (SO->NodeList[3*n1]   + SO->NodeList[3*n2]   + SO->NodeList[3*n3]  )/3; /* centroid of triangle */
         tc[1] = (SO->NodeList[3*n1+1] + SO->NodeList[3*n2+1] + SO->NodeList[3*n3+1])/3; 
         tc[2] = (SO->NodeList[3*n1+2] + SO->NodeList[3*n2+2] + SO->NodeList[3*n3+2])/3;
         norm[0] = SO->FaceNormList[3*i]; norm[1] = SO->FaceNormList[3*i+1]; norm[2] = SO->FaceNormList[3*i+2];
         /* normal is expected to be normalized...*/
         norm[0] *= NormScale; norm[1] *= NormScale; norm[2] *= NormScale;
         fprintf (fout,"%f %f %f \t%f %f %f\n", tc[0], tc[1], tc[2], tc[0]+norm[0], tc[1]+norm[1], tc[2]+norm[2]);            
      }
      
      fclose(fout); fout = NULL;
   }
   
   if (Do_cord) {
      sprintf(OutName, "%s.coord.1D.dset", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Edge output file exists.\nWill not overwrite.");
         exit(1);
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#Cartesian coords, \n");
      fprintf (fout,"#  Center is: [%f %f %f] \n", SO->Center[0], SO->Center[1], SO->Center[2]);
      fprintf (fout,"#nI = Node Index\n");
      fprintf (fout,"#x  = X \n");
      fprintf (fout,"#y  = Y\n");
      fprintf (fout,"#z  = Z\n");
      fprintf (fout,"#nI\tx\ty\tz\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i < SO->N_Node; ++i) {
            fprintf (fout,"%d\t%f\t%f\t%f\t\n",
                  i, SO->NodeList[3*i], SO->NodeList[3*i+1],SO->NodeList[3*i+2]);
      }
      
      fclose(fout); fout = NULL;
   }  
   
   if (Do_edges) {
      
      SUMA_S_Note("Writing edges...");
      
      if (!SO->EL) {
         SUMA_S_Err("Edge list not computed.");
         exit(1);
      }
      
      sprintf(OutName, "%s.edges", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Edge output file exists.\nWill not overwrite.");
         exit(1);
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#Edge List\n");
      fprintf (fout,"#eI = Edge Index\n");
      fprintf (fout,"#n1 = Node 1\n");
      fprintf (fout,"#n2 = Node 2\n");
      fprintf (fout,"#nt = Number of triangles containing edge\n"); 
      fprintf (fout,"#eL = Edge Length\n");
      fprintf (fout,"#eI\tn1\tn2\tnt\teL\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i < SO->EL->N_EL; ++i) {
         if (SO->EL->ELps[i][2] >= 0) {
            n1 = SO->EL->EL[i][0];
            n2 = SO->EL->EL[i][1];
            nt = SO->EL->ELps[i][2];
            n1_3 = 3 * n1;
            n2_3 = 3 * n2;
            edgeL2 = ( (SO->NodeList[n2_3] - SO->NodeList[n1_3]) * (SO->NodeList[n2_3] - SO->NodeList[n1_3]) ) +
                     ( (SO->NodeList[n2_3+1] - SO->NodeList[n1_3+1]) * (SO->NodeList[n2_3+1] - SO->NodeList[n1_3+1]) ) +
                     ( (SO->NodeList[n2_3+2] - SO->NodeList[n1_3+2]) * (SO->NodeList[n2_3+2] - SO->NodeList[n1_3+2]) ); 
                     
            fprintf (fout,"%d\t%d\t%d\t%d\t%f\n",
                  i, n1, n2, nt, sqrt(edgeL2));
                  
         }   
      }
      fclose(fout); fout = NULL;
      
   }
   
   if (Do_area) {
      SUMA_S_Note("Writing areas...");
      
      if (!SO->PolyArea) {
         SUMA_S_Err("Areas not computed");
         exit(1);
      }  
      
      sprintf(OutName, "%s.area", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Area output file exists.\nWill not overwrite.");
         exit(1);
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }
      
      fprintf (fout,"#FaceSet Area\n");
      fprintf (fout,"#fI = FaceSet Index\n");
      fprintf (fout,"#fA = FaceSet Area\n");
      fprintf (fout,"#fI\t#fA\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      for (i=0; i < SO->N_FaceSet; ++i) {
         fprintf (fout,"%d\t%f\n", i, SO->PolyArea[i]);
      }  
      
      fclose(fout); fout = NULL;
   }
   
   if (Do_curv) {
      SUMA_S_Note("Writing curvatures ...");
      
      if (!SO->SC) {
         SUMA_S_Err("Curvatures not computed");
         exit(1);
      }
      
      sprintf(OutName, "%s.curv.1D.dset", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Curvature output file exists.\nWill not overwrite.");
         exit(1);
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }  
      
      fprintf (fout,"#Curvature\n");
      fprintf (fout,"#nI = Node Index\n");
      fprintf (fout,"#T1 = 1 x 3 vector of 1st principal direction of surface\n");
      fprintf (fout,"#T2 = 1 x 3 vector of 2nd principal direction of surface\n");
      fprintf (fout,"#Kp1 = curvature along T1\n");
      fprintf (fout,"#Kp2 = curvature along T2\n");
      fprintf (fout,"#nI\tT1[0]\tT1[1]\tT1[2]\tT2[0]\tT2[1]\tT2[2]\tKp1\tKp2\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i < SO->N_Node; ++i) {
         fprintf (fout,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
            i, SO->SC->T1[i][0], SO->SC->T1[i][1], SO->SC->T1[i][2], 
            SO->SC->T2[i][0], SO->SC->T2[i][1], SO->SC->T2[i][2],
            SO->SC->Kp1[i], SO->SC->Kp2[i] );
      }
      
      fclose(fout); fout = NULL;
   }
   
   if (Do_conv) {
      SUMA_S_Note("Writing convexities ...");
      Cx = (float *)SUMA_GetCx(SO->idcode_str, SUMAg_CF->DsetList, 0);
      if (!Cx) {
         SUMA_S_Err("Convexities not computed");
         exit(1);
      }
      
      sprintf(OutName, "%s.conv.1D.dset", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Convexities output file exists.\nWill not overwrite.");
         exit(1);
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }  
      
      fprintf (fout,"#Convexity\n");
      fprintf (fout,"#nI = Node Index\n");
      fprintf (fout,"#C = Convexity\n");
      fprintf (fout,"#nI\tC\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      for (i=0; i < SO->N_Node; ++i) {
         fprintf (fout,"%d\t%f\n", i, Cx[i]);
      }
      
      fclose(fout); fout = NULL;
   }   
   
   if (Do_vol) {
      float vol;
      fprintf (SUMA_STDOUT,"Calculating surface volume...\n");
      vol = SUMA_Mesh_Volume(SO, NULL, -1);
      fprintf (SUMA_STDERR,   "Volume of closed surface is %f (units3).\n"
                              "Signed volume is  %f (units3).\n", fabs(vol), vol); 
   }
   
   if (Do_en || Do_in) {
      byte *enmask = NULL;
      if (Do_en) fprintf (SUMA_STDOUT,"finding boundary nodes...\n");
      else fprintf (SUMA_STDOUT,"finding internal nodes...\n");
      enmask = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte));  
      if (!enmask) {
         SUMA_S_Crit("Failed to allocate.");
         exit(1);
      }
      if (!SO->EL) {
         SUMA_S_Err("Unexpected error, no edge list.");
         exit(1);
      }
      /* first find the boundary nodes */
      i=0;
      while(i<SO->EL->N_EL) {
         /* find edges that form boundaries */
         if (SO->EL->ELps[i][2] == 1) {
            enmask[SO->EL->EL[i][0]] = 1;
            enmask[SO->EL->EL[i][1]] = 1; 
         }
         ++i;
      }
      
      if (Do_en) {
         sprintf(OutName, "%s.boundarynodes.1D.dset", OutPrefix);
         if (SUMA_filexists(OutName)) {
            SUMA_S_Err("Boundarynodes output file exists.\nWill not overwrite.");
            exit(1);
         }
      } else {
         /* actually need internal ones */
         i=0;
         while(i<SO->EL->N_EL) {
            /* find edges that are internal, NOT part of boundary already */
            if (SO->EL->ELps[i][2] > 1) {
               if (!enmask[SO->EL->EL[i][0]]) enmask[SO->EL->EL[i][0]] = 2;
               if (!enmask[SO->EL->EL[i][1]]) enmask[SO->EL->EL[i][1]] = 2; 
            }
            ++i;
         }
         sprintf(OutName, "%s.internalnodes.1D.dset", OutPrefix);
         if (SUMA_filexists(OutName)) {
            SUMA_S_Err("Internalnodes output file exists.\nWill not overwrite.");
            exit(1);
         }
      }
      
      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }  
      
      if (Do_en) {
         fprintf (fout,"#Boundary Nodes\n");
      } else {
         fprintf (fout,"#Internal Nodes\n");
      }
      fprintf (fout,"#nI = Node Index\n");
      fprintf (fout,"#nI\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);
      
      if (Do_en) {
         for (i=0; i < SO->N_Node; ++i) {
            if (enmask[i] == 1) fprintf (fout,"%d\n", i);
         }
      } else {
         for (i=0; i < SO->N_Node; ++i) {
            if (enmask[i] == 2) fprintf (fout,"%d\n", i);
         }
      }
      SUMA_free(enmask); enmask = NULL;
      
      fclose(fout); fout = NULL;
   }
   
   if (closest_to_xyz) { /* read xyz list, report closest node (SLOW implementation...)*/
      MRI_IMAGE *im = NULL;
      float *far=NULL, *p=NULL;
      int nx2, N_XYZ = 0,i, i3, *closest=NULL, n;
      double xyz[3], d, *dXYZ=NULL;

      /* load the 1D file */
      im = mri_read_1D (closest_to_xyz);
      if (!im) {
         fprintf(SUMA_STDERR,"Error %s:\n Failed to read/find %s.\n", FuncName, closest_to_xyz);
         exit(1);
      }   

      far = MRI_FLOAT_PTR(im);
      if (im->nx == 0) {
         fprintf(SUMA_STDERR,"Error %s:\n Empty file %s.\n", FuncName, closest_to_xyz);
         exit(1);
      }
      if (im->ny != 3) {
         fprintf(SUMA_STDERR,"Error %s:\n Found %d columns in %s. Expecting 3\n", FuncName, im->ny, closest_to_xyz);
         exit(1);
      }

      /* set the results vector */
      dXYZ = (double *)SUMA_malloc(im->nx*sizeof(double));
      closest = (int *)SUMA_malloc(im->nx*sizeof(int));
      if (!dXYZ || !closest) {
         SUMA_S_Crit("Failed to allocate.");
         exit(1);
      }

      nx2 = 2*im->nx;
      for (i=0;i<im->nx; ++i) {  /* for each of the coordinates in question */
         xyz[0] = (double)far[i];
         xyz[1] = (double)far[i+im->nx];
         xyz[2] = (double)far[i+nx2];
         dXYZ[i] = 1023734552736672366372.0;
         closest[i] = -1;
         for (n=0; n<SO->N_Node; ++n) {
            p = &(SO->NodeList[SO->NodeDim*n]);
            SUMA_SEG_NORM(p, xyz, d);
            if (d < dXYZ[i]) {
               dXYZ[i] = d; closest[i] = n;
            }
         }
      }

      /* write out the results */
      sprintf(OutName, "%s.closest.1D.dset", OutPrefix);
      if (SUMA_filexists(OutName)) {
         SUMA_S_Err("Closest nodes output file exists.\nWill not overwrite.");
         exit(1);
      }

      fout = fopen(OutName,"w");
      if (!fout) {
         SUMA_S_Err("Failed to open file for writing.\nCheck your permissions.\n");
         exit(1);
      }  

      fprintf (fout,"#closest nodes reference points\n");
      fprintf (fout,"#n = index of closest node\n");
      fprintf (fout,"#d = distance of closest node\n");
      fprintf (fout,"#X = X of reference point\n");
      fprintf (fout,"#Y = Y of reference point\n");
      fprintf (fout,"#Z = Z of reference point\n");
      fprintf (fout,"#Xn = X of node\n");
      fprintf (fout,"#Yn = Y of node\n");
      fprintf (fout,"#Zn = Z of node\n");
     
      fprintf (fout,"#n\td\tX\tY\tZ\tXn\tYn\tZn\n\n");
      if (histnote) fprintf (fout,"#History:%s\n", histnote);

      for (i=0; i < im->nx; ++i) {
         fprintf (fout,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", 
            closest[i], dXYZ[i], far[i], far[i+im->nx], far[i+nx2], 
               SO->NodeList[SO->NodeDim*closest[i]], 
               SO->NodeList[SO->NodeDim*closest[i]+1], 
               SO->NodeList[SO->NodeDim*closest[i]+2]);
      }

      fclose(fout); fout = NULL;

      /* clean up im */
      if (im) mri_free(im); im = NULL; 

      /* clean up other */
      if (closest) SUMA_free(closest); closest = NULL;
      if (dXYZ) SUMA_free(dXYZ); dXYZ = NULL;
   }
   
   SUMA_LH("Clean up");
   /* clean up */
   if (!SUMA_FreeSpecFields(&Spec)) { SUMA_S_Err("Error freeing"); exit(1); }

   if (MetricList) SUMA_free(MetricList);
   if (OutPrefix) SUMA_free(OutPrefix);
   if (OutName) SUMA_free(OutName);   
   if (SO) SUMA_Free_Surface_Object(SO);
   
   /* dset and its contents are freed in SUMA_Free_CommonFields */
   if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);

   if (histnote) SUMA_free(histnote);
   
   SUMA_RETURN(0);
} /* Main */