/*USE This sample to start writing standalone programs. Change SurfSmooth to the program name of your choosing. */ #include "SUMA_suma.h" SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/ SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< 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 = NULL; /*%< Global pointer to Displayable Object structure vector*/ int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */ SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */ void usage_SUMA_SurfSmooth () { static char FuncName[]={"usage_SUMA_SurfSmooth"}; char * s = NULL, *st = NULL, *sm = NULL; s = SUMA_help_basics(); st = SUMA_help_talk(); sm = SUMA_help_mask(); printf ("\nUsage: SurfSmooth <-spec SpecFile> <-surf_A insurf> <-met method> \n" "\n" " Some methods require additional options detailed below.\n" " I recommend using the -talk_suma option to watch the \n" " progression of the smoothing in real-time in suma.\n" "\n" " Method specific options:\n" " HEAT: <-input inData.1D> <-fwhm F> <-Niter N> \n" " This method is used to filter data\n" " on the surface. It is faster and more stable than\n" " the older LB_FEM below.\n"/* " LB_FEM: <-input mnData.1D> <-fwhm f>\n" " This method is used to filter data\n" " on the surface.\n" */ " LM: [-kpb k] [-lm l m] [-surf_out surfname] [-iw weights]\n" " This method is used to filter the surface's\n" " geometry (node coordinates).\n" " NN_geom: smooth by averaging coordinates of \n" " nearest neighbors.\n" " This method causes shrinkage of surface\n" " and is meant for test purposes only.\n" "\n" " Common options:\n" " [-Niter N] [-output out.1D] [-h/-help] [-dbg_n node]\n" " [-add_index] [-ni_text|-ni_binary] [-talk_suma]\n\n" "\n" " Detailed usage:\n" " -spec SpecFile: Name of specfile containing wurface of interest.\n" " If the surface does not hava 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 now obsolete.\n" " -met method: name of smoothing method to use. Choose from:\n" " HEAT: The newer method by Chung et al. [Ref. 3&4 below]\n" " This method is used for filtering \n" " data on the surface and not for smoothing \n" " the surface's geometry per se. \n" " LB_FEM: (Olde, better use HEAT)\n" " The method by Chung et al. 03.\n" " This method is used for filtering \n" $ " data on the surface and not for smoothing the\n" " surface's geometry per se. See References below.\n" " LM: The smoothing method proposed by G. Taubin 2000\n" " This method is used for smoothing\n" " a surface's geometry. See References below.\n" " NN_geom: A simple nearest neighbor coordinate smoothing.\n" " This interpolation method causes surface shrinkage\n" " that might need to be corrected with the -match_*\n" " options below. \n" "\n" /* " Options for LB_FEM:\n" " It is now recommended that you use the newer method HEAT (see below).\n" " -input inData: file containing data (in 1D or niml format)\n" " $ Each column in inData is processed separately.\n" " The number of rows must equal the number of\n" " nodes in the surface. You$can select certain\n" " columns using the [] notation adopted by AFNI's\n" " programs.\n" " Note: The program will infer the format of the input\n" " file from the extension of inData. \n" " -fwhm f: Full Width at Half Maximum in surface coordinate units (usuallly mm)\n" " of an equivalent Gaussian filter had the surface been flat.\n" " With curved surfaces, the equation used to estimate FWHM is \n" " an approximation. \n" " Blurring on the surface depends on the geodesic instead \n" " of the Euclidean distances. See Ref #1 for more details \n" " on this parameter.\n" */ "\n" " Options for HEAT:\n" " -input inData : file containing data (in 1D or NIML format)\n" " Each column in inData is processed separately.\n" " The number of rows must equal the number of\n" " nodes in the surface. You can select certain\n" " columns using the [] notation adopted by AFNI's\n" " Note: The program will infer the format of the input\n" " file from the extension of inData. \n" " programs.\n" " -fwhm F: Effective Full Width at Half Maximum in surface coordinate units \n" " (usuallly mm) of an equivalent Gaussian filter had the surface been \n" " flat. With curved surfaces, the equation used to estimate FWHM is \n" " an approximation. For Gaussian filters, FWHM, SIGMA (STD-DEV) and RMS\n" " FWHM = 2.354820 * SIGMA = 1.359556 * RMS\n" " Blurring on the surface depends on the geodesic instead \n" " of the Euclidean distances. \n" " and one of the following two parameters:\n" " (See Refs #3&4 for more details)\n" " -Niter N: Number of iterations (default is -1).\n" " You can now set this parameter to -1 and have \n" " the program suggest a value based on the -fwhm value.\n" " Too large or too small a number of iterations can affect \n" " smoothing results. Acceptable values depend on \n" " the average distance between nodes oj the mesh and\n" " the desired fwhm.\n" " -sigma S: Bandwidth of smoothing kernel (for a single iteration).\n" " S shoul` be small (< 1) and is related to the previous two\n" " parameters by: F = sqrt(N) * S * 2.355\n" "\n" "\n" " Options for LM:\n" " -kpb k: Band pass frequency (default is 0.1).\n" " values should be in the range 0 < k < 10\n" " -lm and -kpb options are mutually exclusive.\n" " -lm l m: Lambda and Mu parameters. Sample values are:\n" " 0.6307 and -.6732\n" " NOTE: -lm and -kpb options are mutually exclusive.\n" " -surf_out surfname: Writes the surface with smoothed coordinates\n" " to disk. For SureFit and 1D formats, only the\n" " coord file is written out.\n" " NOTE: -surf_out and -output are mutually exclusive.\n" " -iw wgt: Set interpolation weights to wgt. You can choose from:\n" " Equal : Equal weighting, fastest (default), \n" " tends to make edges equal.\n" " Fujiwara: Weighting based on inverse edge length.\n" " Would be a better preserver of geometry when\n" " mesh has irregular edge lengths.\n" " Desbrun : Weighting based on edge angles (slooow).\n" " Removes tangential displacement during smoothing.\n" " Might not be too useful for brain surfaces.\n" "\n" " Options for NN_geom:\n" " -match_size r: Adjust node coordinates of smoothed surface to \n" " approximates the original's size.\n" " Node i on the filtered surface is repositioned such$\n" " that |c i| = 1/N sum(|cr j|) where\n" " c and cr are the centers of the smoothed and original\n" " surfaces, respectively.\n" " N is the number of nodes that are within r [surface \n" " coordinate units] along the surface (geodesic) from node i.\n" " j is one of the nodes neighboring i.\n" " -match_vol tol: Adjust node coordinates of smoothed surface to \n" " approximates the original's volume.\n" " Nodes on the filtered surface are repositioned such\n" " that the volume of the filtered surface equals, \n" " withmn tolerence tol, that of the original surface. \n" " See option -vol in SurfaceMetrics for information about\n" " and calculation of the volume of a closed surface.\n" " -match_area tol: Adjust node coordinates of smoothed surface to \n" " approximates the original's surface.\n" " Nodes on the filtered surface are repositioned such\n" " that the surface of the filtered surface equals, \n" " within tolerance tol, that of the original surface. \n" " -match_sphere rad: Project nodes of smoothed surface to a sphere\n" " of radius rad. Projection is carried out along the \n" " direction formed by the surface's center and the node.\n" " -surf_out surfname: Writes the surface with smoothed coordinates\n" " to disk. For SureFit and 1D formats, only the\n" " coord file is written out.\n" "\n" " Common options:\n" " -Niter N: Number of smoothing iterations (default is 100)\n" " For practical reasons, this number must be a multiple of 2\n" " NOTE 1: For HEAT method, you can set Niter to -1, in conjunction\n" " with -fwhm FWHM option, and the program\n" " will pick an acceptable number for you.\n" " NOTE 2: For LB_FEM method, the number of iterations controls the\n" " iteration steps (dt in Ref #1).\n" " dt = fwhm*fwhm / (16*Niter*log(2));\n" " dt must satisfy conditions that depend on the internodal\n" " distance and the spatial derivatives of the signals being \n" " filtered on the surface.\n" " As a rule of thumb, if increasing Niter does not alter\n" " the results then your choice is fine (smoothing has converged).\n" " For an example of the artifact caused by small Niter see:\n" " http://afni.nimh.nih.gov/sscc/staff/ziad/SUMA/SuSmArt/DSart.html\n" " To avoid this problem altogether, it is better that you use \n" " the newer method HEAT which does not suffer from this problem.\n" " -output OUT: Name of output file. \n" " The default is inData_sm with LB_FEM and HEAT method\n" " and NodeList_sm with LM method.\n" " NOTE: For data smoothing methods like HEAT, If a format extension,\n" " such as .1D.dset or .niml.dset is present in OUT, \n" " then the output will be written in that format. \n" " Otherwise, the format is the same as the input's\n" " -add_index : Output the node index in the first column.\n" " This is not done by default.\n" " -dbg_n node : output debug information for node 'node'.\n" " -use_neighbors_outside_mask: When using -c_mask or -b_mask or -n_mask\n" " options, allow value from a node nj neighboring\n" " node n to contribute to the value at n even if\n" " nj is not in the mask.\n" " The default is to ignore all nodes not in the\n" " mask.\n" "\n" "%s" "\n" "%s" "\n" "%s" "\n" " Sample commands lines for using SurfSmooth:\n" " The surface used in this example had no spec file, so \n" " a quick.spec was created using:\n" " quickspec -tn 1D NodeList.1D FaceSetList.1D \n" "\n" " Sample commands lines for data smoothing:\n" " SurfSmooth -spec quick.spec -surf_A NodeList.1D -met HEAT \\\n" " -input in.1D -Niter -1 -fwhm 8 -add_index \\\n" " -output in_smh8.1D.dset \n" /* " Or using the older (less recommended method):\n" " SurfSmooth -spec quick.spec -surf_A NodeList.1D -met LB_FEM \\\n" " -input in.1D -Niter 100 -fwhm 8 -add_index \\\n" " -output in_sm8.1D.dset \n" " This command filters (on the surface) the data in in.1D\n" " and puts the output in in_sm8.1D.dset with the first column \n" " containing the node index and the second containing the \n" " filtered version of in.1D.\n" " \n"*/ "\n" " You can colorize the input and output data using ScaleToMap:\n" " ScaleToMap -input in.1D 0 1 -cmap BGYR19 \\\n" " -clp MIN MAX > in.1D.col \\\n" " ScaleToMap -input in_sm8.1D 0 1 -cmap BGYR19 \\\n" " -clp MIN MAX > in_sm8.1D.col \\\n" "\n" " For help on using ScaleToMap see ScaleToMap -help\n" " Note that the MIN MAX represent the minimum and maximum\n" " values in in.1D. You should keep them constant in both \n" " commands in order to be able to compare the resultant colorfiles.\n" " You can import the .col files with the 'c' command in SUMA.\n" "\n" " You can send the data to SUMA with each iteration.\n" " To do so, start SUMA with these options:\n" " suma -spec quick.spec -niml &\n" " and add these options to SurfSmooth's command line above:\n" " -talk_suma -refresh_rate 5\n" "\n" " Sample commands lines for surface smoothing:\n" " SurfSmooth -spec quick.spec -surf_A NodeList.1D -met LM \\\n" " -output NodeList_sm100.1D -Niter 100 -kpb 0.1 \n" " This command smoothes the surface's geometry. The smoothed\n" " node coordinates are written out to NodeList_sm100.1D. \n" "\n" " Sample command for considerable surface smoothing and inflation\n" " back to original volume:\n" " SurfSmooth -spec quick.spec -surf_A NodeList.1D -met NN_geom \\\n" " -output NodeList_inflated_mvol.1D -Niter 1500 \\\n" " -match_vol 0.01\n" " Sample command for considerable surface smoothing and inflation\n" " back to original area:\n" " SurfSmooth -spec quick.spec -surf_A NodeList.1D -met NN_geom \\\n" " -output NodeList_inflated_marea.1D -Niter 1500 \\\n" " -match_area 0.01\n" "\n" " References: \n" " (1) M.K. Chung et al. Deformation-based surface morphometry\n" " applied to gray matter deformation. \n" " Neuroimage 18 (2003) 198-213\n" " M.K. Chung Statistical morphometry in computational\n" " neuroanatomy. Ph.D. thesis, McGill Univ.,\n" " Montreal, Canada\n" " (2) G. Taubin. Mesh Signal Processing. \n" " Eurographics 2000.\n" " (3) M.K. Chung et al. Cortical thickness analysis in autism \n" " via heat kernel smoothing. NeuroImage, \n" " submitted.(2005) \n" " http://www.stat.wisc.edu/~mchung/papers/ni_heatkernel.pdf\n" " (4) M.K. Chung, Heat kernel smoothing and its application to \n" " cortical manifolds. Technical Report 1090. \n" " Department of Statististics, U.W.Madison\n" " http://www.stat.wisc.edu/~mchung/papers/heatkernel_tech.pdf" "\n" " See Also: \n" " ScaleToMap to colorize the output, however it is better\n" " to load surface datasets directly into SUMA and colorize\n" " them interactively." "\n" "\n", sm, st, s); SUMA_free(s); s = NULL; SUMA_free(st); st = NULL; SUMA_free(sm); sm = 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"); exit (0); } #define SURFSMOOTH_MAX_SURF 1 /*!< Maximum number of input surfaces */ typedef enum { SUMA_NO_METH, SUMA_LB_FEM_1D, SUMA_LB_FEM, SUMA_LM, SUMA_BRUTE_FORCE, SUMA_NN_GEOM, SUMA_HEAT_05_1D, SUMA_HEAT_05} SUMA_SMOOTHING_METHODS; typedef struct { float OffsetLim; float lim; float fwhm; float kpb; float l; float m; int ShowNode; int Method; int dbg; int N_iter; int AddIndex; int insurf_method; /* method used to specify input surfaces. 0 then none input 1 the old way 2 the new (-spec way) */ SUMA_SO_File_Type iType; char *vp_name; char *sv_name; char *if_name; char *if_name2; char *in_name; char *out_name; /* this one's dynamically allocated so you'll have to free it yourself */ char *ShowOffset_DBG; /* meant to be a file where one outputs some debugging info. Not being used ...*/ char *surf_out; char *surf_names[SURFSMOOTH_MAX_SURF]; char *spec_file; int MatchMethod; byte *nmask; byte strict_mask; float sigma; SUMA_DSET_FORMAT oform; SUMA_GENERIC_ARGV_PARSE *ps; } SUMA_SURFSMOOTH_OPTIONS; /*! \brief parse the arguments for SurfSmooth program \param argv (char *) \param argc (int) \return Opt (SUMA_SURFSMOOTH_OPTIONS *) options structure. To free it, use SUMA_free(Opt->out_name); SUMA_free(Opt); */ SUMA_SURFSMOOTH_OPTIONS *SUMA_SurfSmooth_ParseInput (char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps) { static char FuncName[]={"SUMA_SurfSmooth_ParseInput"}; SUMA_SURFSMOOTH_OPTIONS *Opt=NULL; int kar, i, ind; char *outname; SUMA_Boolean brk = NOPE; SUMA_Boolean LocalHead = NOPE; SUMA_ENTRY; Opt = (SUMA_SURFSMOOTH_OPTIONS *)SUMA_malloc(sizeof(SUMA_SURFSMOOTH_OPTIONS)); kar = 1; Opt->OffsetLim = -1.0; Opt->MatchMethod = 0; Opt->lim = 1000000.0; Opt->fwhm = -1; Opt->ShowNode = -1; Opt->Method = SUMA_NO_METH; Opt->dbg = 0; Opt->if_name = NULL; Opt->if_name2 = NULL; Opt->in_name = NULL; Opt->out_name = NULL; Opt->vp_name = NULL; Opt->sv_name = NULL; Opt->surf_out = NULL; Opt->ShowOffset_DBG = NULL; Opt->iType = SUMA_FT_NOT_SPECIFIED; Opt->N_iter = -100; Opt->kpb = -1.0; Opt->l = -1.0; Opt->m = -1.0; Opt->AddIndex = 0; Opt->insurf_method = 0; Opt->nmask = NULL; Opt->spec_file = NULL; Opt->sigma = -1.0; Opt->oform = SUMA_NO_DSET_FORMAT; Opt->strict_mask = 1; Opt->ps = ps; SUMA_Set_Taubin_Weights(SUMA_EQUAL); for (i=0; isurf_names[i] = NULL; } outname = NULL; brk = NOPE; 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_SurfSmooth(); exit (0); } SUMA_SKIP_COMMON_OPTIONS(brk, kar); if (!brk && strcmp(argv[kar], "-dist") == 0) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -dist \n"); exit (1); } Opt->OffsetLim = atof(argv[kar]); if (Opt->OffsetLim <= 0 && Opt->OffsetLim != -1.0) { fprintf (SUMA_STDERR, "Bad value (%f) for refresh_rate\n", Opt->OffsetLim); exit (1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-dbg_n") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need 1 argument after -dbg_n \n"); exit (1); } Opt->ShowNode = atoi(argv[kar]); SUMA_Set_SurfSmooth_NodeDebug(Opt->ShowNode); brk = YUP; } if (!brk && (strcmp(argv[kar], "-Niter") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need 1 arguments after -Niter \n"); exit (1); } Opt->N_iter = atoi(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-kpb") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need 1 arguments after -kpb \n"); exit (1); } if (Opt->l != -1.0 || Opt->m != -1.0) { fprintf (SUMA_STDERR, "options -lm and -kpb are mutually exclusive\n"); exit (1); } Opt->kpb = atof(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-surf_out") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need 1 arguments after -surf_out\n"); exit (1); } if (outname) { fprintf (SUMA_STDERR, "-output and -surf_out are mutually exclusive.\n"); exit(1); } Opt->surf_out = argv[kar]; brk = YUP; } if (!brk && (strcmp(argv[kar], "-lm") == 0)) { kar ++; if (kar+1 >= argc) { fprintf (SUMA_STDERR, "need 2 arguments after -lm \n"); exit (1); } if (Opt->kpb != -1.0) { fprintf (SUMA_STDERR, "options -lm and -kpb are mutually exclusive\n"); exit (1); } Opt->l = atof(argv[kar]); kar ++; Opt->m = atof(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-iw") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -iw \n"); exit (1); } SUMA_TO_LOWER(argv[kar]); if (strcmp(argv[kar], "equal") == 0 ) { SUMA_Set_Taubin_Weights(SUMA_EQUAL); }else if (strcmp(argv[kar],"fujiwara")==0) { SUMA_Set_Taubin_Weights(SUMA_FUJIWARA); }else if (strcmp(argv[kar],"desbrun")==0) { SUMA_Set_Taubin_Weights(SUMA_DESBRUN); } else { fprintf (SUMA_STDERR, "Weights option %s not understood.\n", argv[kar]); exit (1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-input") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -input \n"); exit (1); } Opt->in_name = argv[kar]; brk = YUP; } if (!brk && (strcmp(argv[kar], "-output") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -output\n"); exit (1); } if (Opt->surf_out) { fprintf (SUMA_STDERR, "options -surf_out and -output are mutually exclusive\n"); exit (1); } outname = argv[kar]; brk = YUP; } if (!brk && (strcmp(argv[kar], "-add_index") == 0)) { Opt->AddIndex = 1; brk = YUP; } if (!brk && (strcmp(argv[kar], "-use_neighbors_outside_mask") == 0)) { Opt->strict_mask = 0; brk = YUP; } if (!brk && (strcmp(argv[kar], "-i_fs") == 0)) { SUMA_SL_Err("Option -i_fs is obsolete.\nUse -spec and -surf_A instead.\n"); exit(1); kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -i_fs \n"); exit (1); } Opt->if_name = argv[kar]; Opt->iType = SUMA_FREE_SURFER; if (!Opt->insurf_method) Opt->insurf_method = 1; else { fprintf (SUMA_STDERR, "already specified input surface.\n"); exit(1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-i_sf") == 0)) { SUMA_SL_Err("Option -i_sf is obsolete.\nUse -spec and -surf_A instead.\n"); exit(1); kar ++; if (kar+1 >= argc) { fprintf (SUMA_STDERR, "need 2 argument after -i_sf\n"); exit (1); } Opt->if_name = argv[kar]; kar ++; Opt->if_name2 = argv[kar]; Opt->iType = SUMA_SUREFIT; if (!Opt->insurf_method) Opt->insurf_method = 1; else { fprintf (SUMA_STDERR, "already specified input surface.\n"); exit(1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-i_vec") == 0)) { SUMA_SL_Err("Option -i_vec is obsolete.\nUse -spec and -surf_A instead.\n"); exit(1); kar ++; if (kar+1 >= argc) { fprintf (SUMA_STDERR, "need 2 argument after -i_vec\n"); exit (1); } Opt->if_name = argv[kar]; kar ++; Opt->if_name2 = argv[kar]; Opt->iType = SUMA_VEC; if (!Opt->insurf_method) Opt->insurf_method = 1; else { fprintf (SUMA_STDERR, "already specified input surface.\n"); exit(1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-i_dx") == 0)) { SUMA_SL_Err("Option -i_ply is obsolete.\nUse -spec and -surf_A instead.\n"); exit(1); kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -i_dx\n "); exit (1); } Opt->if_name = argv[kar]; Opt->iType = SUMA_OPENDX_MESH; if (!Opt->insurf_method) Opt->insurf_method = 1; else { fprintf (SUMA_STDERR, "already specified input surface.\n"); exit(1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-i_ply") == 0)) { SUMA_SL_Err("Option -i_ply is obsolete.\nUse -spec and -surf_A instead.\n"); exit(1); kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -i_ply\n "); exit (1); } Opt->if_name = argv[kar]; Opt->iType = SUMA_PLY; if (!Opt->insurf_method) Opt->insurf_method = 1; else { fprintf (SUMA_STDERR, "already specified input surface.\n"); exit(1); } brk = YUP; } if (!brk && (strcmp(argv[kar], "-spec") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -spec \n"); exit (1); } Opt->spec_file = argv[kar]; if (!Opt->insurf_method) Opt->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 >= SURFSMOOTH_MAX_SURF) { fprintf (SUMA_STDERR, "-surf_X SURF_NAME option is out of range,\n" " only surf_A allowed.\n"); exit (1); } kar ++; Opt->surf_names[ind] = argv[kar]; if (Opt->insurf_method && Opt->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], "-lim") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -lim \n"); exit (1); } Opt->lim = atof(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-match_size") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -match_size \n"); exit (1); } Opt->lim = atof(argv[kar]); Opt->MatchMethod = 1; brk = YUP; } if (!brk && (strcmp(argv[kar], "-match_vol") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -match_vol \n"); exit (1); } Opt->lim = atof(argv[kar]); Opt->MatchMethod = 2; brk = YUP; } if (!brk && (strcmp(argv[kar], "-match_area") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -match_area \n"); exit (1); } Opt->lim = atof(argv[kar]); Opt->MatchMethod = 3; brk = YUP; } if (!brk && (strcmp(argv[kar], "-match_sphere") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -match_sphere \n"); exit (1); } Opt->lim = atof(argv[kar]); Opt->MatchMethod = 4; brk = YUP; } if (!brk && (strcmp(argv[kar], "-fwhm") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -fwhm \n"); exit (1); } Opt->fwhm = atof(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-sigma") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -sigma \n"); exit (1); } Opt->sigma = atof(argv[kar]); brk = YUP; } if (!brk && (strcmp(argv[kar], "-met") == 0)) { kar ++; if (kar >= argc) { fprintf (SUMA_STDERR, "need argument after -met \n"); exit (1); } if (strcmp(argv[kar], "LB_FEM_1D") == 0) Opt->Method = SUMA_LB_FEM_1D; else if (strcmp(argv[kar], "LB_FEM") == 0) Opt->Method = SUMA_LB_FEM; else if (strcmp(argv[kar], "LM") == 0) Opt->Method = SUMA_LM; else if (strcmp(argv[kar], "BF") == 0) Opt->Method = SUMA_BRUTE_FORCE; else if (strcmp(argv[kar], "NN_geom") == 0) Opt->Method = SUMA_NN_GEOM; else if (strcmp(argv[kar], "HEAT_1D") == 0) Opt->Method = SUMA_HEAT_05_1D; else if (strcmp(argv[kar], "HEAT") == 0) Opt->Method = SUMA_HEAT_05; else { fprintf (SUMA_STDERR, "Method %s not supported.\n", argv[kar]); exit (1); } if (Opt->Method == SUMA_HEAT_05) { if (Opt->N_iter < 0) { Opt->N_iter = -1; } } else { if (Opt->N_iter < 0) { Opt->N_iter = -Opt->N_iter; } } brk = YUP; } if (!brk && !ps->arg_checked[kar]) { fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]); exit (1); } else { brk = NOPE; kar ++; } } /* check on options for HEAT budiness first */ if (Opt->Method == SUMA_HEAT_05_1D || Opt->Method == SUMA_HEAT_05) { float sequiv; if (Opt->N_iter < 0 && Opt->fwhm > 0 && Opt->sigma < 0) { /* will make a suggestion, but after surface is read */ }else { if ( (Opt->N_iter < 0 && Opt->fwhm < 0 && Opt->sigma < 0) || (Opt->N_iter > 0 && Opt->fwhm > 0 && Opt->sigma > 0) || (Opt->N_iter < 0 && Opt->fwhm < 0) || (Opt->N_iter < 0 && Opt->sigma < 0)|| (Opt->fwhm < 0 && Opt->sigma < 0) ) { fprintf (SUMA_STDERR,"Error %s:\n" "Need to specify two and only two out of the three options:\n" "-N_iter, -fwhm, -sigma\n", FuncName); exit (1); } if (Opt->N_iter < 0) { sequiv = Opt->fwhm * 0.42466090; Opt->N_iter = SUMA_POW2(sequiv/Opt->sigma); if (Opt->N_iter % 2) ++Opt->N_iter; if (Opt->N_iter < 100) { fprintf (SUMA_STDERR,"Warning: N_iter = %d\n" "Perhaps too small for comfort.\n" "Consider reducing sigma.\n", Opt->N_iter); } /* recalculate sigma */ sequiv = Opt->fwhm * 0.42466090; Opt->sigma = sequiv / sqrt(Opt->N_iter); } if (Opt->fwhm < 0) { sequiv = sqrt(Opt->N_iter)*Opt->sigma; Opt->fwhm = sequiv / 0.42466090; } if (Opt->sigma < 0) { sequiv = Opt->fwhm * 0.42466090; Opt->sigma = sequiv / sqrt(Opt->N_iter); } } } if (Opt->Method != SUMA_HEAT_05) { if (Opt->N_iter == -1) { /* default */ Opt->N_iter = 100; } if (Opt->N_iter < 1) { fprintf (SUMA_STDERR,"Error %s:\nWith -Niter N option, N must be > 1\n", FuncName); exit (1); } } else { if (Opt->N_iter == -1) {Opt->N_iter = -2; }/* let it slide, we'll take care of it after surface is loaded in main */ } if ( (Opt->N_iter % 2) && (Opt->Method == SUMA_LB_FEM_1D || Opt->Method == SUMA_LB_FEM || Opt->Method == SUMA_HEAT_05_1D || Opt->Method == SUMA_LM) ) { fprintf (SUMA_STDERR, "Number of iterations must be a multiple of 2.\n%d is not a multiple of 2.\n", Opt->N_iter); exit(1); } if (Opt->ShowNode < 0 && Opt->ShowOffset_DBG) { fprintf (SUMA_STDERR,"Error %s:\nBad debug node index (%d) in option -dbg_n\n", FuncName, Opt->ShowNode); exit (1); } if (Opt->insurf_method == 1) { SUMA_SL_Err("Obsolete method for surface specification.\nShould not have gotten here."); exit(1); } /* can't test for file existence here because of square brackets */ if (0 && Opt->in_name && !SUMA_filexists(Opt->in_name)) { fprintf (SUMA_STDERR,"Error %s:\n%s not found.\n", FuncName, Opt->if_name); exit(1); } if (Opt->Method == SUMA_NO_METH) { fprintf (SUMA_STDERR,"Error %s:\nNo method was specified.\n", FuncName); exit(1); } if (ps->cs->talk_suma && Opt->insurf_method != 2) { fprintf (SUMA_STDERR, "must specify surface using -spec option\n" "if you whish to talk to suma.\n"); exit(1); } if (Opt->insurf_method == 2) { if (!Opt->surf_names[0] || !Opt->spec_file) { fprintf (SUMA_STDERR, "failed to specify either -spec or -surf_X options.\n"); exit(1); } } if (outname) { if (SUMA_filexists(outname)) { fprintf (SUMA_STDERR,"Error %s:\noutput file %s exists.\n", FuncName, outname); exit(1); } Opt->out_name = SUMA_copy_string(outname); Opt->oform = SUMA_GuessFormatFromExtension(Opt->out_name); } else { switch (Opt->Method) { case SUMA_LB_FEM_1D: /* form autoname */ Opt->out_name = SUMA_Extension(Opt->in_name, ".1D", YUP); /*remove .1D */ Opt->out_name = SUMA_append_replace_string(Opt->out_name,"_sm", "", 1); /* add _sm to prefix */ Opt->out_name = SUMA_append_replace_string(Opt->out_name,".1D", "", 1); /* add .1D */ break; case SUMA_LB_FEM: /* form autoname */ Opt->oform = SUMA_GuessFormatFromExtension(Opt->in_name); Opt->out_name = SUMA_RemoveDsetExtension_s(Opt->in_name, Opt->oform); Opt->out_name = SUMA_append_replace_string(Opt->out_name, "_sm", "", 1); /* add _sm to prefix */ Opt->out_name = SUMA_append_replace_string(Opt->out_name, (char*)SUMA_ExtensionOfDsetFormat (Opt->oform), "", 1); /* add extension */ break; case SUMA_LM: /* form autoname */ Opt->out_name = SUMA_copy_string("NodeList_sm.1D"); break; case SUMA_BRUTE_FORCE: /* form autoname */ Opt->out_name = SUMA_copy_string("NodeList_Offsetsm.1D"); break; case SUMA_NN_GEOM: /* form autoname */ Opt->out_name = SUMA_copy_string("NodeList_NNsm.1D"); break; case SUMA_HEAT_05_1D: /* form autoname */ Opt->out_name = SUMA_Extension(Opt->in_name, ".1D", YUP); /*remove .1D */ Opt->out_name = SUMA_append_replace_string(Opt->out_name,"_smh", "", 1); /* add _smh to prefix */ Opt->out_name = SUMA_append_replace_string(Opt->out_name,".1D", "", 1); /* add .1D */ break; case SUMA_HEAT_05: /* form autoname */ Opt->oform = SUMA_GuessFormatFromExtension(Opt->in_name); Opt->out_name = SUMA_RemoveDsetExtension_s(Opt->in_name, Opt->oform); Opt->out_name = SUMA_append_replace_string(Opt->out_name,"_sm", "", 1); /* add _sm to prefix */ Opt->out_name = SUMA_append_replace_string(Opt->out_name, (char*)SUMA_ExtensionOfDsetFormat (Opt->oform), "", 1); /* add extension */ break; default: fprintf (SUMA_STDERR,"Error %s:\nNot ready for this option here.\n", FuncName); exit(1); break; } if (SUMA_filexists(Opt->out_name)) { fprintf (SUMA_STDERR,"Error %s:\noutput file %s exists.\n", FuncName, Opt->out_name); exit(1); } } /* method specific checks */ switch (Opt->Method) { case SUMA_LB_FEM_1D: case SUMA_LB_FEM: if (!Opt->in_name) { fprintf (SUMA_STDERR,"Error %s:\ninput data not specified.\n", FuncName); exit(1); } if (Opt->fwhm == -1.0) { fprintf (SUMA_STDERR,"Error %s:\n-fwhm option must be used with -met LB_FEM.\n", FuncName); exit(1); }else if (Opt->fwhm <= 0.0) { fprintf (SUMA_STDERR,"Error %s:\nFWHM must be > 0\n", FuncName); exit(1); } if (Opt->kpb >= 0) { fprintf (SUMA_STDERR,"Error %s:\n-kpb option is not valid with -met LB_FEM.\n", FuncName); exit(1); } break; case SUMA_HEAT_05: case SUMA_HEAT_05_1D: if (!Opt->in_name) { fprintf (SUMA_STDERR,"Error %s:\ninput data not specified.\n", FuncName); exit(1); } if (Opt->fwhm == -1.0) { fprintf (SUMA_STDERR,"Error %s:\n-fwhm option must be used with -met HEAT.\n", FuncName); exit(1); }else if (Opt->fwhm <= 0.0) { fprintf (SUMA_STDERR,"Error %s:\nFWHM must be > 0\n", FuncName); exit(1); } if (Opt->kpb >= 0) { fprintf (SUMA_STDERR,"Error %s:\n-kpb option is not valid with -met HEAT.\n", FuncName); exit(1); } break; case SUMA_BRUTE_FORCE: break; case SUMA_LM: if ( (Opt->l != -1.0 || Opt->m != -1.0) && Opt->kpb != -1.0) { fprintf (SUMA_STDERR,"Error %s:\nYou cannot mix options -kpb and -lm \n", FuncName); exit(1); } if (Opt->kpb != -1.0 && (Opt->kpb < 0.000001 || Opt->kpb > 10)) { fprintf (SUMA_STDERR,"Error %s:\nWith -kpb k option, you should satisfy 0 < k < 10\n", FuncName); exit(1); } if (Opt->l == -1.0 && Opt->m == -1.0 && Opt->kpb == -1.0) { Opt->kpb = 0.1; } if (Opt->l == -1.0 || Opt->m == -1.0) { /* convert kpb into l and m */ if (!SUMA_Taubin_Smooth_Coef (Opt->kpb, &(Opt->l), &(Opt->m))) { SUMA_SL_Err("Failed to find smoothing coefficients"); exit(1); } } if (Opt->in_name) { fprintf (SUMA_STDERR,"Error %s:\nOption -input not valid with -met LM.\n", FuncName); exit(1); } if (Opt->fwhm != -1.0) { fprintf (SUMA_STDERR,"Error %s:\nOption -fwhm not valid with -met LM.\n", FuncName); exit(1); } break; case SUMA_NN_GEOM: if (Opt->in_name) { fprintf (SUMA_STDERR,"Error %s:\nOption -input not valid with -met NN_geom.\n", FuncName); exit(1); } if (0 && Opt->lim > 1000) { fprintf (SUMA_STDERR,"Error %s:\n-lim option not specified.\n", FuncName); exit(1); } if (Opt->fwhm != -1.0) { fprintf (SUMA_STDERR,"Error %s:\nOption -fwhm not valid with -met NN_geom.\n", FuncName); exit(1); } break; default: fprintf (SUMA_STDERR,"Error %s:\nNot ready for this option here.\n", FuncName); exit(1); break; } if (0 && ((Opt->ps->bmaskname && Opt->ps->nmaskname) || (Opt->ps->bmaskname && Opt->ps->cmask) || (Opt->ps->nmaskname && Opt->ps->cmask) ) ) { fprintf (SUMA_STDERR,"Error %s:\n-n_mask, -b_mask, and -c_mask options are mutually exclusive.\n", FuncName); exit(1); }else { /* SUMA_S_Warn("For testing! Turn me back on!!!\n"); */ /* Now it is allowed */ } SUMA_RETURN (Opt); } int main (int argc,char *argv[]) {/* Main */ static char FuncName[]={"SurfSmooth"}; int kar, icol, nvec, ncol=0, i, ii, N_inmask = -1; float *data_old = NULL, *far = NULL; float **DistFirstNeighb; void *SO_name = NULL; SUMA_SurfaceObject *SO = NULL, *SOnew = NULL; MRI_IMAGE *im = NULL; SUMA_SFname *SF_name = NULL; struct timeval start_time, start_time_all; float etime_GetOffset, etime_GetOffset_all; SUMA_GET_OFFSET_STRUCT *OffS = NULL; SUMA_SURFSMOOTH_OPTIONS *Opt; FILE *fileout=NULL; float **wgt=NULL, *dsmooth=NULL; SUMA_INDEXING_ORDER d_order=SUMA_NO_ORDER; SUMA_COMM_STRUCT *cs = NULL; SUMA_SurfSpecFile Spec; SUMA_GENERIC_ARGV_PARSE *ps=NULL; SUMA_DSET *dset = NULL; int iform; SUMA_Boolean LocalHead = NOPE; SUMA_STANDALONE_INIT; SUMA_mainENTRY; /* Allocate space for DO structure */ SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS); ps = SUMA_Parse_IO_Args(argc, argv, "-o;-i;-sv;-talk;-mask;"); if (argc < 6) { usage_SUMA_SurfSmooth(); exit (1); } Opt = SUMA_SurfSmooth_ParseInput (argv, argc, ps); cs = ps->cs; if (!cs) exit(1); /* now for the real work */ if (Opt->insurf_method == 1) { /* method 1 */ SUMA_SL_Err("Input in this method is no longer supported.\n"); exit(1); } else { /* method 2 */ int SO_read = -1; if (!SUMA_AllocSpecFields(&Spec)) { SUMA_S_Err("Failed to allocate"); exit(1); } if (!SUMA_Read_SpecFile (Opt->spec_file, &Spec)) { fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName); exit(1); } SO_read = SUMA_spec_select_surfs(&Spec, Opt->surf_names, SURFSMOOTH_MAX_SURF, 0); if ( SO_read != 1 ) { fprintf(SUMA_STDERR,"Error %s: Found %d surfaces, expected only 1.\n", FuncName, SO_read); exit(1); } /* now read into SUMAg_DOv */ if (!SUMA_LoadSpec_eng(&Spec, SUMAg_DOv, &SUMAg_N_DOv, Opt->sv_name, 0, SUMAg_CF->DsetList) ) { fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_LoadSpec_eng\n", FuncName); exit(1); } /* now identify surface needed */ SO = SUMA_find_named_SOp_inDOv(Opt->surf_names[0], SUMAg_DOv, SUMAg_N_DOv); } if (!SO) { fprintf (SUMA_STDERR,"Error %s: Failed to read input surface.\n", FuncName); exit (1); } /* setup the mask, if needed */ if (Opt->nmask) { SUMA_S_Err("Should be null here!"); exit(1); } if (!(Opt->nmask = SUMA_load_all_command_masks(Opt->ps->bmaskname, Opt->ps->nmaskname, Opt->ps->cmask, SO->N_Node, &N_inmask)) && N_inmask < 0) { SUMA_S_Err("Failed loading mask"); exit(1); } if (Opt->nmask) { fprintf(SUMA_STDOUT,"%d nodes in mask:\n", N_inmask); if (LocalHead) { ii = 0; for (i=0; iN_Node; ++i) { if (Opt->nmask[i]) { fprintf(SUMA_STDERR,"%6d ", i); ++ii; if (!(ii % 12)) fprintf(SUMA_STDERR,"\n"); } } fprintf(SUMA_STDERR,"\n"); } } else { if (LocalHead) { fprintf(SUMA_STDOUT,"No masking.\n"); } } if (Opt->ShowNode >= 0 && Opt->ShowNode >= SO->N_Node) { fprintf (SUMA_STDERR,"Error %s: Requesting debugging info for a node index (%d) \n" "that does not exist in a surface of %d nodes.\nRemember, indexing starts at 0.\n", FuncName, Opt->ShowNode, SO->N_Node); exit (1); } /* form EL and FN */ if (!SO->EL || !SO->FN) { /* normally you'd call SUMA_SurfaceMetrics_eng (SO, "EdgeList", NULL, 0) but that should be done in SUMA_SurfaceMetrics_eng. */ SUMA_SLP_Err("Unexpexted NULL SO->EL or SO->FN"); exit(1); } /* see if SUMA talk is turned on */ if (ps->cs->talk_suma) { cs->istream = SUMA_GEOMCOMP_LINE; if ( Opt->Method == SUMA_LB_FEM_1D || Opt->Method == SUMA_LB_FEM || Opt->Method == SUMA_HEAT_05_1D || Opt->Method == SUMA_HEAT_05 ) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NO_DSET_TYPE, 0)) { SUMA_SL_Err("Failed to initialize SUMA_SendToSuma"); cs->Send = NOPE; ps->cs->talk_suma = NOPE; } }else if (Opt->Method == SUMA_LM) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NO_DSET_TYPE, 0)) { SUMA_SL_Err("Failed to initialize SUMA_SendToSuma"); cs->Send = NOPE; ps->cs->talk_suma = NOPE; } }else if (Opt->Method == SUMA_NN_GEOM) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NO_DSET_TYPE, 0)) { SUMA_SL_Err("Failed to initialize SUMA_SendToSuma"); cs->Send = NOPE; ps->cs->talk_suma = NOPE; } }else if (Opt->Method == SUMA_BRUTE_FORCE) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NO_DSET_TYPE, 0)) { SUMA_SL_Err("Failed to initialize SUMA_SendToSuma"); cs->Send = NOPE; ps->cs->talk_suma = NOPE; } }else { SUMA_SL_Err("Can't talk to suma with the chosen method.\n"); ps->cs->talk_suma = NOPE; } } ncol = 3; /* default for geometry smoothing. That is changed below for data smoothing. */ switch (Opt->Method) { case SUMA_HEAT_05_1D:/* Operates on 1D files, OBSOLETE but still accessible with -met HEAT_1D */ /* Moo Chung's method for interpolation weights */ { if (Opt->N_iter < 0) { SUMA_S_Errv("Bad number of iterations (%d)\n", Opt->N_iter); exit(1); } /* now load the input data */ im = mri_read_1D (Opt->in_name); if (!im) { SUMA_SL_Err("Failed to read 1D file"); exit(1); } far = MRI_FLOAT_PTR(im); nvec = im->nx; ncol = im->ny; d_order = SUMA_COLUMN_MAJOR; if (!nvec) { SUMA_SL_Err("Empty file"); exit(1); } if (nvec != SO->N_Node) { fprintf(SUMA_STDERR, "Error %s:\n" "Expecting 1D file to have %d rows\n" " found %d rows instead.\n", FuncName, SO->N_Node, nvec); exit(1); } if (LocalHead) SUMA_etime(&start_time,0); wgt = SUMA_Chung_Smooth_Weights_05(SO, Opt->sigma); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: weight computation took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } dsmooth = SUMA_Chung_Smooth_05 ( SO, wgt, Opt->N_iter, Opt->fwhm, far, ncol, SUMA_COLUMN_MAJOR, NULL, cs, Opt->nmask, Opt->strict_mask); if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL; } break; case SUMA_LB_FEM_1D: /* Operates on 1D files, OBSOLETE but still accessible with -met LB_FEM_1D */ /* Moo Chung's method for interpolation weights */ { /* now load the input data */ im = mri_read_1D (Opt->in_name); if (!im) { SUMA_SL_Err("Failed to read 1D file"); exit(1); } far = MRI_FLOAT_PTR(im); nvec = im->nx; ncol = im->ny; d_order = SUMA_COLUMN_MAJOR; if (!nvec) { SUMA_SL_Err("Empty file"); exit(1); } if (nvec != SO->N_Node) { fprintf(SUMA_STDERR, "Error %s:\n" "Expecting 1D file to have %d rows\n" " found %d rows instead.\n", FuncName, SO->N_Node, nvec); exit(1); } if (LocalHead) SUMA_etime(&start_time,0); wgt = SUMA_Chung_Smooth_Weights(SO); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: weight computation took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } dsmooth = SUMA_Chung_Smooth ( SO, wgt, Opt->N_iter, Opt->fwhm, far, ncol, SUMA_COLUMN_MAJOR, NULL, cs, Opt->nmask, Opt->strict_mask); if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } #if 0 /* writing is now done below ... */ fileout = fopen(Opt->out_name, "w"); if (Opt->AddIndex) SUMA_disp_vecmat (dsmooth, SO->N_Node, ncol, 1, d_order, fileout, YUP); else SUMA_disp_vecmat (dsmooth, SO->N_Node, ncol, 1, d_order, fileout, NOPE); fclose(fileout); fileout = NULL; if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL; #endif if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL; } break; case SUMA_LB_FEM: /* Moo Chung's method for interpolation weights with dsets */ { /* now load the input data */ iform = SUMA_GuessFormatFromExtension(Opt->in_name); if (!(dset = SUMA_LoadDset_s (Opt->in_name, &iform, 0))) { SUMA_S_Err("Failed to read dset"); exit(1); } if (Opt->oform == SUMA_NO_DSET_FORMAT) Opt->oform = iform; if (!SDSET_VECLEN(dset) || !SDSET_VECNUM(dset)) { SUMA_SL_Err("Empty file"); exit(1); } if (SDSET_VECLEN(dset) != SO->N_Node) { if (LocalHead) fprintf(SUMA_STDERR, "Warning %s:\n" "Expecting 1D file to have %d rows\n" " found %d rows instead.\n" "Function should deal with this properly but check results\n", FuncName, SO->N_Node, SDSET_VECLEN(dset)); } if (Opt->AddIndex || Opt->oform == SUMA_NIML) { if (!SUMA_AddNodeIndexColumn(dset, SO->N_Node)) { SUMA_S_Err("Failed to add a node index column"); exit(1); } } if (LocalHead) SUMA_etime(&start_time,0); wgt = SUMA_Chung_Smooth_Weights(SO); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: weight computation took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (!SUMA_Chung_Smooth_dset ( SO, wgt, Opt->N_iter, Opt->fwhm, dset, cs, Opt->nmask, Opt->strict_mask)) { SUMA_S_Err("Failed in SUMA_Chung_Smooth_dset"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL; } break; case SUMA_HEAT_05: /* Moo Chung's method for interpolation weights */ { double avg_wt, sequiv; if (Opt->N_iter < 0 && Opt->fwhm > 0 && Opt->sigma < 0) { /* make a suggestion */ Opt->sigma = sqrt(-SO->EL->AvgLe/(2*log(0.01))); /* making the average SUMA_CHUNG_KERNEL_NUMER be 1 percent */ /* have sigma and fwhm, what is N? */ sequiv = Opt->fwhm * 0.42466090; Opt->N_iter = SUMA_POW2(sequiv/Opt->sigma); if (Opt->N_iter % 2) ++Opt->N_iter; if (Opt->N_iter < 4) Opt->N_iter = 4; /* need a few iterations */ /* now reset sigma based on number of iterations */ sequiv = Opt->fwhm * 0.42466090; Opt->sigma = sequiv / sqrt(Opt->N_iter); } else if (Opt->N_iter < 0) { SUMA_S_Errv("Negative number of iterations (%d).Should not be here.\n", Opt->N_iter); exit(1); } fprintf (SUMA_STDERR,"Effective FWHM = %f, kernel bandwidth = %f, N_iter = %d\n", Opt->fwhm, Opt->sigma, Opt->N_iter); /* check that sigma is not too small or too big relative to average segment length */ avg_wt = SUMA_CHUNG_KERNEL_NUMER(SO->EL->AvgLe,Opt->sigma); fprintf(SUMA_STDERR, "Kernel Bandwidth / Average Edge Distance = %f/%f = %f\n" " Corresponding Kernel Numerator = %g\n", Opt->sigma, SO->EL->AvgLe, Opt->sigma/SO->EL->AvgLe, avg_wt); if (avg_wt > 0.05) { SUMA_S_Warnv("Average weight assigned per node is of %g\n" "It is advisable to increase the number of\n" "iterations, if possible.\n", avg_wt); } if (avg_wt < 1e-4) { SUMA_S_Warnv("Average weight assigned per node is of %g\n" "It is advisable to decrease the number of\n" "iterations, if possible.\n", avg_wt); } if (avg_wt < 1e-7) { SUMA_S_Warnv("Average weight assigned per node is of %g\n" "Cannot trust results. Reduce number of\n" "iterations, if possible.\n", avg_wt); } /* now load the input data */ iform = SUMA_GuessFormatFromExtension(Opt->in_name); if (!(dset = SUMA_LoadDset_s (Opt->in_name, &iform, 0))) { SUMA_S_Err("Failed to read dset"); exit(1); } if (LocalHead) { SUMA_LHv("Input %s:\n", Opt->in_name); SUMA_ShowDset(dset, 0, NULL); } if (Opt->oform == SUMA_NO_DSET_FORMAT) Opt->oform = iform; if (!SDSET_VECLEN(dset) || !SDSET_VECNUM(dset)) { SUMA_SL_Err("Empty file"); exit(1); } if (SDSET_VECLEN(dset) != SO->N_Node) { if (LocalHead) fprintf(SUMA_STDERR, "Warning %s:\n" "Expecting 1D file to have %d rows\n" " found %d rows instead.\n" "Function should deal with this properly but check results\n", FuncName, SO->N_Node, SDSET_VECLEN(dset)); } if (Opt->AddIndex || Opt->oform == SUMA_NIML) { if (!SUMA_AddNodeIndexColumn(dset, SO->N_Node)) { SUMA_S_Err("Failed to add a node index column"); exit(1); } } if (LocalHead) SUMA_etime(&start_time,0); wgt = SUMA_Chung_Smooth_Weights_05(SO, Opt->sigma); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: weight computation took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (!SUMA_Chung_Smooth_05_dset ( SO, wgt, Opt->N_iter, Opt->fwhm, dset, cs, Opt->nmask, Opt->strict_mask)) { SUMA_S_Err("Failed in SUMA_Chung_Smooth_05_dset"); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL; } break; case SUMA_BRUTE_FORCE: { /* now load the input data */ iform = SUMA_GuessFormatFromExtension(Opt->in_name); if (!(dset = SUMA_LoadDset_s (Opt->in_name, &iform, 0))) { SUMA_S_Err("Failed to read dset"); exit(1); } if (Opt->oform == SUMA_NO_DSET_FORMAT) Opt->oform = iform; if (!SDSET_VECLEN(dset) || !SDSET_VECNUM(dset)) { SUMA_SL_Err("Empty file"); exit(1); } if (SDSET_VECLEN(dset) != SO->N_Node) { if (LocalHead) fprintf(SUMA_STDERR, "Warning %s:\n" "Expecting 1D file to have %d rows\n" " found %d rows instead.\n" "Function should deal with this properly but check results\n", FuncName, SO->N_Node, SDSET_VECLEN(dset)); } if (Opt->AddIndex || Opt->oform == SUMA_NIML) { if (!SUMA_AddNodeIndexColumn(dset, SO->N_Node)) { SUMA_S_Err("Failed to add a node index column"); exit(1); } } if (LocalHead) SUMA_etime(&start_time,0); if (LocalHead) { SUMA_LHv("Input %s:\n", Opt->in_name); SUMA_ShowDset(dset, 0, NULL); } if (!SUMA_Offset_Smooth_dset( SO, Opt->fwhm, Opt->OffsetLim, Opt->N_iter, dset, cs, Opt->nmask, Opt->strict_mask)) { SUMA_S_Err("Failed in SUMA_Offset_Smooth_dset "); exit(1); } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (0 && LocalHead) { SUMA_LH("See dsmooth.1D"); SUMA_disp_vecmat (dsmooth, SO->N_Node, 3, 1, d_order, NULL, YUP); } /* writing of results is done below */ etime_GetOffset_all = SUMA_etime(&start_time_all,1); fprintf(SUMA_STDERR, "%s: Done.\nSearch to %f mm took %f minutes for %d nodes.\n" , FuncName, Opt->lim, etime_GetOffset_all / 60.0 , SO->N_Node); } break; case SUMA_NN_GEOM: /* brute forcem nearset neighbor interpolation */ { if (LocalHead) SUMA_etime(&start_time,0); d_order = SUMA_ROW_MAJOR; dsmooth = SUMA_NN_GeomSmooth( SO, Opt->N_iter, SO->NodeList, 3, d_order, NULL, cs, Opt->nmask, Opt->strict_mask); if (0 && LocalHead) { SUMA_LH("See dsmooth.1D"); SUMA_disp_vecmat (dsmooth, SO->N_Node, 3, 1, d_order, NULL, YUP); } if (!dsmooth) { SUMA_SL_Err("Failed in SUMA_NN_Geom_Smooth"); exit(1); } /* writing of results is done below */ } break; case SUMA_LM: /* Taubin's */ { if (LocalHead) SUMA_etime(&start_time,0); wgt = NULL; d_order = SUMA_ROW_MAJOR; if (SUMA_Get_Taubin_Weights() == SUMA_FUJIWARA) { SUMA_SL_Note("Fujiwara!!!"); wgt = SUMA_Taubin_Fujiwara_Smooth_Weights(SO, NULL, NULL); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } } else if (SUMA_Get_Taubin_Weights() == SUMA_DESBRUN) { wgt = SUMA_Taubin_Desbrun_Smooth_Weights(SO, NULL, NULL); if (!wgt) { SUMA_SL_Err("Failed to compute weights.\n"); exit(1); } } else if (SUMA_Get_Taubin_Weights() != SUMA_EQUAL) { /* fprintf(stderr,"%d, %d\n", SUMA_Get_Taubin_Weights() , SUMA_EQUAL); */ SUMA_SL_Err("Weights improperly initialized!"); exit(1); } dsmooth = SUMA_Taubin_Smooth (SO, wgt, Opt->l, Opt->m, SO->NodeList, Opt->N_iter, 3, d_order, NULL, cs, Opt->nmask, Opt->strict_mask); if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL; /* writing of results is done below */ } break; default: SUMA_SL_Err("Bad method, should not be here."); exit(1); break; } if (Opt->Method == SUMA_NN_GEOM) { if (Opt->MatchMethod) { SUMA_LH("Fixing shrinkage..."); /* Create a surface that is a descendant of SO, use the new coordinates */ SOnew = SUMA_CreateChildSO( SO, dsmooth, SO->N_Node, NULL, -1, 0); /* SOnew->NodeList is now == dsmooth */ /* fix the shrinking ..*/ switch (Opt->MatchMethod) { case 1: if (!SUMA_EquateSurfaceSize(SOnew, SO, Opt->lim, cs)) { SUMA_SL_Warn("Failed to fix surface size.\nTrying to finish ..."); } /* send the unshrunk bunk */ if (cs->Send) { SUMA_LH("Sending last fix to SUMA ..."); if (!SUMA_SendToSuma (SO, cs, (void *)SOnew->NodeList, SUMA_NODE_XYZ, 1)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted."); } } /* to make matters parallel with the other methods, keep dsmooth and free SOnew */ dsmooth = SOnew->NodeList; /* CHECK IF THAT's the case here... coordinates have a new pointer after Equating surface size */ SOnew->NodeList = NULL; /* new coordinates will stay alive in dsmooth */ SUMA_Free_Surface_Object(SOnew); SOnew=NULL; break; case 2: if (!SUMA_EquateSurfaceVolumes(SOnew, SO, Opt->lim, cs)) { SUMA_SL_Warn("Failed to fix surface size.\nTrying to finish ..."); } /* send the unshrunk bunk */ if (cs->Send) { SUMA_LH("Sending last fix to SUMA ..."); if (!SUMA_SendToSuma (SO, cs, (void *)SOnew->NodeList, SUMA_NODE_XYZ, 1)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted."); } } /* to make matters parallel with the other methods, keep dsmooth and free SOnew */ dsmooth = SOnew->NodeList; /* coordinates have a new pointer after Equating surface volumes */ SOnew->NodeList = NULL; /* new coordinates will stay alive in dsmooth */ SUMA_Free_Surface_Object(SOnew); SOnew=NULL; break; case 3: if (!SUMA_EquateSurfaceAreas(SOnew, SO, Opt->lim, cs)) { SUMA_SL_Warn("Failed to fix surface size.\nTrying to finish ..."); } /* send the unshrunk bunk */ if (cs->Send) { SUMA_LH("Sending last fix to SUMA ..."); if (!SUMA_SendToSuma (SO, cs, (void *)SOnew->NodeList, SUMA_NODE_XYZ, 1)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted."); } } /* to make matters parallel with the other methods, keep dsmooth and free SOnew */ dsmooth = SOnew->NodeList; /* coordinates have a new pointer after Equating surface volumes */ SOnew->NodeList = NULL; /* new coordinates will stay alive in dsmooth */ SUMA_Free_Surface_Object(SOnew); SOnew=NULL; break; case 4: if (!SUMA_ProjectSurfaceToSphere(SOnew, SO, Opt->lim, cs)) { SUMA_SL_Warn("Failed to fix surface size.\nTrying to finish ..."); } /* send the unshrunk bunk */ if (cs->Send) { SUMA_LH("Sending last fix to SUMA ..."); if (!SUMA_SendToSuma (SO, cs, (void *)SOnew->NodeList, SUMA_NODE_XYZ, 1)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted."); } } /* to make matters parallel with the other methods, keep dsmooth and free SOnew */ dsmooth = SOnew->NodeList; /* CHECK IF THAT's the case here... coordinates have a new pointer after Equating surface size */ SOnew->NodeList = NULL; /* new coordinates will stay alive in dsmooth */ SUMA_Free_Surface_Object(SOnew); SOnew=NULL; break; case 0: break; default: SUMA_SL_Err("Huh ?"); break; } } if (LocalHead) { etime_GetOffset = SUMA_etime(&start_time,1); fprintf(SUMA_STDERR, "%s: Total processing took %f seconds for %d nodes.\n" "Projected time per 100000 nodes is: %f minutes\n", FuncName, etime_GetOffset, SO->N_Node, etime_GetOffset * 100000 / 60.0 / (SO->N_Node)); } } /* write out the filtered geometry. Should not be executed for data smoothing */ if (Opt->surf_out) { if (!dsmooth) { SUMA_SL_Err("NULL dsmooth for geometry smoothing. Either failed to smooth or logical error."); exit(1); } SUMA_free(SO->NodeList); SO->NodeList = dsmooth; dsmooth = NULL; /* replace NodeList */ switch (SO->FileType) { case SUMA_SUREFIT: SF_name = (SUMA_SFname *) SUMA_malloc(sizeof(SUMA_SFname)); sprintf(SF_name->name_coord,"%s", Opt->surf_out); SF_name->name_topo[0] = '\0'; SO_name = (void *)SF_name; if (!SUMA_Save_Surface_Object (SO_name, SO, SUMA_SUREFIT, SUMA_ASCII, NULL)) { fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName); exit (1); } break; case SUMA_VEC: SF_name = (SUMA_SFname *) SUMA_malloc(sizeof(SUMA_SFname)); sprintf(SF_name->name_coord,"%s", Opt->surf_out); SF_name->name_topo[0] = '\0'; SO_name = (void *)SF_name; if (!SUMA_Save_Surface_Object (SO_name, SO, SUMA_VEC, SUMA_ASCII, NULL)) { fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName); exit (1); } break; case SUMA_FREE_SURFER: SO_name = (void *)Opt->surf_out; if (!SUMA_Save_Surface_Object (SO_name, SO, SUMA_FREE_SURFER, SUMA_ASCII, NULL)) { fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName); exit (1); } break; case SUMA_FREE_SURFER_PATCH: fprintf (SUMA_STDERR,"Error %s: No support for writing Free Surfer patches.\n", FuncName); exit (1); break; case SUMA_PLY: SO_name = (void *)Opt->surf_out; if (!SUMA_Save_Surface_Object (SO_name, SO, SUMA_PLY, SUMA_FF_NOT_SPECIFIED, NULL)) { fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName); exit (1); } break; case SUMA_OPENDX_MESH: SO_name = (void *)Opt->surf_out; if (!SUMA_Save_Surface_Object (SO_name, SO, SUMA_OPENDX_MESH, SUMA_ASCII, NULL)) { fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName); exit (1); } break; default: fprintf (SUMA_STDERR,"Error %s: Bad format.\n", FuncName); exit(1); } } else { if (Opt->Method != SUMA_LB_FEM && Opt->Method != SUMA_HEAT_05 && Opt->Method != SUMA_BRUTE_FORCE) { if (!dsmooth) { SUMA_SL_Err("NULL dsmooth for data smoothing. Either failed to smooth or logical error."); exit(1); } fileout = fopen(Opt->out_name, "w"); if (Opt->AddIndex) SUMA_disp_vecmat (dsmooth, SO->N_Node, ncol, 1, d_order, fileout, YUP); else SUMA_disp_vecmat (dsmooth, SO->N_Node, ncol, 1, d_order, fileout, NOPE); fclose(fileout); fileout = NULL; } else if (Opt->Method == SUMA_LB_FEM || Opt->Method == SUMA_HEAT_05 || Opt->Method == SUMA_BRUTE_FORCE) { SUMA_NEWDSET_ID_LABEL_HIST(dset, Opt->out_name) ; if (Opt->AddIndex) SUMA_SetAddIndex_1D(1); SUMA_WriteDset_s(Opt->out_name, dset, Opt->oform, 0, 0); SUMA_FreeDset(dset); dset = NULL; } else { SUMA_S_Err("Fix me"); } } /* you don't want to exit rapidly because the SUMA might not be done processing the last elements*/ if (cs->Send && !cs->GoneBad) { /* cleanup and close connections */ if (Opt->Method == SUMA_LB_FEM_1D) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NODE_RGBAb, 2)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCleanup failed"); } }else if (Opt->Method == SUMA_LM) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NODE_XYZ, 2)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCleanup failed"); } }else if (Opt->Method == SUMA_NN_GEOM) { if (!SUMA_SendToSuma (SO, cs, NULL, SUMA_NODE_XYZ, 2)) { SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCleanup failed"); } } } SUMA_LH("clean up"); if (SO->NodeList != dsmooth) { SUMA_LH("Freeing dsmooth...:"); if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL; SUMA_LH("Done:"); } if (!SUMA_FreeSpecFields(&Spec)) { SUMA_S_Err("Failed to free spec fields"); } mri_free(im); im = NULL; /* done with that baby */ if (cs) cs = NULL; /* ps->cs if freed below */ if (SF_name) SUMA_free(SF_name); if (Opt->insurf_method == 1) { if (SO) SUMA_Free_Surface_Object(SO); } if (data_old) SUMA_free(data_old); if (Opt->out_name) SUMA_free(Opt->out_name); Opt->out_name = NULL; if (Opt->nmask) SUMA_free(Opt->nmask); Opt->nmask = NULL; if (Opt) SUMA_free(Opt); if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL; if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) { SUMA_SL_Err("DO Cleanup Failed!"); } exit(0); }