Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

SUMA_SurfToSurf.c

Go to the documentation of this file.
00001 #include "SUMA_suma.h"
00002 
00003 SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/
00004 SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00005                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00006 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00007 SUMA_DO *SUMAg_DOv = NULL;   /*!< Global pointer to Displayable Object structure vector*/
00008 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00009 SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */
00010 
00011 void usage_SurfToSurf (SUMA_GENERIC_ARGV_PARSE *ps)
00012 {
00013       static char FuncName[]={"usage_SurfToSurf"};
00014       char * s = NULL, *sio=NULL, *st = NULL, *sts = NULL;
00015       int i;
00016       s = SUMA_help_basics();
00017       sio  = SUMA_help_IO_Args(ps);
00018       printf ( "\n"
00019                "Usage: SurfToSurf <-i_TYPE S1> [<-sv SV1>]\n"
00020                "                  <-i_TYPE S2> [<-sv SV1>]\n"
00021                "                  [<-prefix PREFIX>]\n"
00022                "                  [<-output_params PARAM_LIST>]\n"
00023                "                  [<-node_indices NODE_INDICES>]\n"
00024                "                  [<-proj_dir PROJ_DIR>]\n"
00025                "                  [<-data DATA>]\n"
00026                "                  [<-node_debug NODE>]\n"
00027                "                  [<-debug DBG_LEVEL>]\n"
00028                "                  [-make_consistent]\n"
00029                " \n"
00030                "  This program is used to interpolate data from one surface (S2)\n"
00031                " to another (S1), assuming the surfaces are quite similar in\n"
00032                " shape but having different meshes (non-isotopic).\n"
00033                " This is done by projecting each node (nj) of S1 along the normal\n"
00034                " at nj and finding the closest triangle t of S2 that is intersected\n"
00035                " by this projection. Projection is actually bidirectional.\n"
00036                " If such a triangle t is found, the nodes (of S2) forming it are \n"
00037                " considered to be the neighbors of nj.\n" 
00038                " Values (arbitrary data, or coordinates) at these neighboring nodes\n"
00039                " are then transfered to nj using barycentric interpolation or \n"
00040                " nearest-node interpolation.\n"
00041                " Nodes whose projections fail to intersect triangles in S2 are given\n"
00042                " nonsensical values of -1 and 0.0 in the output.\n" 
00043                "\n"
00044                " Mandatory input:\n"
00045                "  Two surfaces are required at input. See -i_TYPE options\n"
00046                "  below for more information. \n"
00047                "\n"
00048                " Optional input:\n"
00049                "  -prefix PREFIX: Specify the prefix of the output file.\n"
00050                "                  The output file is in 1D format at the moment.\n"
00051                "                  Default is SurfToSurf\n"
00052                "  -output_params PARAM_LIST: Specify the list of mapping\n"
00053                "                             parameters to include in output\n"
00054                "     PARAM_LIST can have any or all of the following:\n"
00055                "        NearestTriangleNodes: Use Barycentric interpolation (default)\n"
00056                "                              and output indices of 3 nodes from S2\n"
00057                "                              that neighbor nj of S1\n"
00058                "        NearestNode: Use only the closest node from S2 (of the three \n"
00059                "                     closest neighbors) to nj of S1 for interpolation\n"
00060                "                     and output the index of that closest node.\n"
00061                "        NearestTriangle: Output index of triangle t from S2 that\n"
00062                "                         is the closest to nj along its projection\n"
00063                "                         direction. \n"  
00064                "        DistanceToSurf: Output distance (signed) from nj, along \n"
00065                "                        projection direction to S2.\n"
00066                "                        This is the parameter output by the precursor\n"
00067                "                        program CompareSurfaces\n"
00068                "        ProjectionOnSurf: Output coordinates of projection of nj onto \n"
00069                "                          triangle t of S2.\n"
00070                "        Data: Output the data from S2, interpolated onto S1\n"
00071                "              If no data is specified via the -data option, then\n"
00072                "              the XYZ coordinates of SO2's nodes are considered\n"
00073                "              the data.\n"
00074                "  -data DATA: 1D file containing data to be interpolated.\n"
00075                "              Each row i contains data for node i of S2.\n"
00076                "              You must have one row for each node making up S2.\n"
00077                "              In other terms, if S2 has N nodes, you need N rows\n"
00078                "              in DATA. \n"
00079                "              Each column of DATA is processed separately (think\n"
00080                "              sub-bricks, and spatial interpolation).\n"
00081                "              You can use [] selectors to choose a subset \n"
00082                "              of columns.\n"
00083                "              If -data option is not specified and Data is in PARAM_LIST\n"
00084                "              then the XYZ coordinates of SO2's nodes are the data.\n"
00085                "  -node_indices NODE_INDICES: 1D file containing the indices of S1\n"
00086                "                              to consider. The default is all of the\n"
00087                "                              nodes in S1. Only one column of values is\n"
00088                "                              allowed here, use [] selectors to choose\n"
00089                "                              the column of node indices if NODE_INDICES\n"
00090                "                              has multiple columns in it.\n"
00091                "  -proj_dir PROJ_DIR: 1D file containing projection directions to use\n"
00092                "                      instead of the node normals of S1.\n"
00093                "                      Each row should contain one direction for each\n"
00094                "                      of the nodes forming S1.\n"
00095                "  -make_consistent: Force a consistency check and correct triangle \n"
00096                "                    orientation of S1 if needed. Triangles are also\n"
00097                "                    oriented such that the majority of normals point\n"
00098                "                    away from center of surface.\n"
00099                "                    The program might not succeed in repairing some\n"
00100                "                    meshes with inconsistent orientation.\n" 
00101                "\n"
00102                "%s"
00103                "%s"
00104                "\n", sio,  s);
00105       SUMA_free(s); s = NULL; SUMA_free(st); st = NULL; SUMA_free(sio); sio = NULL;       
00106       s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
00107       printf("       Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov     \n"
00108              "       Shruti Japee LBC/NIMH/NIH  shruti@codon.nih.gov \n");
00109       exit(0);
00110 }
00111 
00112 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_SurfToSurf_ParseInput(char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
00113 {
00114    static char FuncName[]={"SUMA_BrainWrap_ParseInput"}; 
00115    SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
00116    int kar;
00117    SUMA_Boolean brk, accepting_out;
00118    
00119    SUMA_Boolean LocalHead = NOPE;
00120 
00121    SUMA_ENTRY;
00122    
00123    Opt = SUMA_Alloc_Generic_Prog_Options_Struct();
00124    kar = 1;
00125    brk = NOPE;
00126    Opt->in_1D = NULL;
00127    Opt->NodeDbg = -1;
00128    Opt->debug = 0;
00129    Opt->NearestNode = 3;
00130    Opt->NearestTriangle = 0;
00131    Opt->DistanceToMesh = 0;
00132    Opt->ProjectionOnMesh = 0;
00133    Opt->Data = 0;
00134    Opt->in_name = NULL;
00135    Opt->out_prefix = NULL;
00136    Opt->fix_winding = 0;
00137    accepting_out = NOPE;
00138    while (kar < argc) { /* loop accross command ine options */
00139                 /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
00140       
00141       if (!brk && accepting_out) { /* make sure you have not begun with new options */
00142          if (*(argv[kar]) == '-') accepting_out = NOPE;
00143       }
00144                 
00145       if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
00146                          usage_SurfToSurf(ps);
00147           exit (0);
00148                 }
00149                 
00150                 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
00151       
00152       if (!brk && (strcmp(argv[kar], "-debug") == 0))
00153       {
00154          if (kar+1 >= argc)
00155          {
00156             fprintf (SUMA_STDERR, "need a number after -debug \n");
00157             exit (1);
00158          }
00159          
00160          Opt->debug = atoi(argv[++kar]);
00161          brk = YUP;
00162       }
00163       
00164       if (!brk && (strcmp(argv[kar], "-node_debug") == 0))
00165       {
00166          if (kar+1 >= argc)
00167          {
00168             fprintf (SUMA_STDERR, "need a number after -node_debug \n");
00169             exit (1);
00170          }
00171          
00172          Opt->NodeDbg = atoi(argv[++kar]);
00173          brk = YUP;
00174       }
00175       
00176       if (!brk && (strcmp(argv[kar], "-node_indices") == 0))
00177       {
00178          if (kar+1 >= argc)
00179          {
00180             fprintf (SUMA_STDERR, "need a parameter after -node_indices \n");
00181             exit (1);
00182          }
00183          
00184          Opt->in_nodeindices = argv[++kar];
00185          brk = YUP;
00186       }
00187       
00188       if (!brk && (strcmp(argv[kar], "-output_params") == 0))
00189       {
00190          if (kar+1 >= argc)
00191          {
00192             fprintf (SUMA_STDERR, "need at least one parameter after output_params \n");
00193             exit (1);
00194          }
00195          
00196          accepting_out = YUP;
00197          brk = YUP;
00198       }
00199       
00200       if (!brk && accepting_out && (strcmp(argv[kar], "NearestNode") == 0)) {
00201          if (Opt->NearestNode < 1) Opt->NearestNode = 1;
00202          brk = YUP;
00203       }
00204       
00205       if (!brk && accepting_out && (strcmp(argv[kar], "NearestTriangleNodes") == 0)) {
00206          if (Opt->NearestNode < 3) Opt->NearestNode = 3;
00207          brk = YUP;
00208       }
00209       
00210       if (!brk && accepting_out && (strcmp(argv[kar], "NearestTriangle") == 0)) {
00211          if (Opt->NearestTriangle < 1) Opt->NearestTriangle = 1;
00212          brk = YUP;
00213       }
00214       
00215       if (!brk && accepting_out && (strcmp(argv[kar], "DistanceToSurf") == 0)) {
00216          Opt->DistanceToMesh = 1;
00217          brk = YUP;
00218       }
00219       
00220       if (!brk && accepting_out && (strcmp(argv[kar], "ProjectionOnSurf") == 0)) {
00221          Opt->ProjectionOnMesh = 1;
00222          brk = YUP;
00223       }
00224       
00225       if (!brk && accepting_out && (strcmp(argv[kar], "Data") == 0)) {
00226          Opt->Data = 1;
00227          brk = YUP;
00228       }
00229       
00230       if (!brk && (strcmp(argv[kar], "-data") == 0))
00231       {
00232          if (kar+1 >= argc)
00233          {
00234             fprintf (SUMA_STDERR, "need a name after -data \n");
00235             exit (1);
00236          }
00237          ++kar;
00238          if (strcmp(argv[kar],"_XYZ_") == 0) { 
00239             /* default Opt->in_name = NULL*/
00240             if (Opt->in_name) {
00241                SUMA_SL_Err("Unexpected non null value");
00242                exit (1);
00243             }
00244          } else {
00245             Opt->in_name = SUMA_copy_string(argv[kar]);
00246          }
00247          Opt->Data = 1;
00248          brk = YUP;
00249       }
00250       
00251       if (!brk && (strcmp(argv[kar], "-prefix") == 0))
00252       {
00253          if (kar+1 >= argc)
00254          {
00255             fprintf (SUMA_STDERR, "need a name after -prefix \n");
00256             exit (1);
00257          }
00258          
00259          Opt->out_prefix = SUMA_Extension(argv[++kar],".1D", YUP);
00260          brk = YUP;
00261       }
00262       
00263       if (!brk && (strcmp(argv[kar], "-proj_dir") == 0))
00264       {
00265          if (kar+1 >= argc)
00266          {
00267             fprintf (SUMA_STDERR, "need a name after -proj_dir \n");
00268             exit (1);
00269          }
00270          
00271          Opt->in_1D = argv[++kar];
00272          brk = YUP;
00273       }
00274       if (!brk && (strcmp(argv[kar], "-make_consistent") == 0))
00275       {
00276          Opt->fix_winding = 1;
00277          brk = YUP;
00278       }
00279       if (!brk && !ps->arg_checked[kar]) {
00280                         fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00281                         exit (1);
00282                 } else {        
00283                         brk = NOPE;
00284                         kar ++;
00285                 }
00286    }
00287    
00288    if (!Opt->out_prefix) Opt->out_prefix = SUMA_copy_string("SurfToSurf"); 
00289    
00290    SUMA_RETURN(Opt);
00291 }
00292 
00293 int main (int argc,char *argv[])
00294 {/* Main */    
00295    static char FuncName[]={"SurfToSurf"}; 
00296    SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;  
00297    SUMA_GENERIC_ARGV_PARSE *ps=NULL;
00298    SUMA_SurfaceObject *SO1=NULL, *SO2 = NULL;
00299    SUMA_SurfSpecFile *Spec = NULL;
00300    SUMA_M2M_STRUCT *M2M = NULL;
00301    int N_Spec=0, *nodeind = NULL, N_nodeind, icol, i, j;
00302    MRI_IMAGE *im = NULL, *im_data=NULL;
00303         int nvec, ncol=0, nvec_data, ncol_data, Nchar;
00304    float *far = NULL, *far_data=NULL, *dt = NULL, *projdir=NULL;
00305    char *outname = NULL, *s=NULL, sbuf[100];
00306    void *SO_name = NULL;   
00307    FILE *outptr=NULL;
00308    SUMA_Boolean exists = NOPE;
00309    SUMA_INDEXING_ORDER d_order;
00310    SUMA_STRING *SS=NULL;
00311    SUMA_Boolean LocalHead = NOPE;
00312 
00313         SUMA_mainENTRY;
00314    SUMA_STANDALONE_INIT;
00315 
00316    /* Allocate space for DO structure */
00317         SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
00318    ps = SUMA_Parse_IO_Args(argc, argv, "-i;-t;-spec;-s;-sv;-o;");
00319    
00320    if (argc < 2) {
00321       usage_SurfToSurf(ps);
00322       exit (1);
00323    }
00324    
00325    Opt = SUMA_SurfToSurf_ParseInput (argv, argc, ps);
00326 
00327    /* if output surface requested, check on pre-existing file */
00328    if (ps->o_N_surfnames) {
00329       SO_name = SUMA_Prefix2SurfaceName(ps->o_surfnames[0], NULL, NULL, ps->o_FT[0], &exists);
00330       if (exists) {
00331          fprintf(SUMA_STDERR,"Error %s:\nOutput file(s) %s* on disk.\nWill not overwrite.\n", FuncName, ps->o_surfnames[0]);
00332          exit(1);
00333       }
00334    } 
00335    
00336    if (Opt->debug > 2) LocalHead = YUP;
00337    
00338    outname = SUMA_append_string(Opt->out_prefix,".1D");
00339    if (SUMA_filexists(outname)) {
00340       fprintf(SUMA_STDERR,"Output file %s exists.\n", outname);
00341       exit(1);
00342    }
00343    
00344    /* Load the surfaces from command line*/
00345    Spec = SUMA_IO_args_2_spec(ps, &N_Spec);
00346    if (N_Spec != 1) {
00347       SUMA_S_Err( "Multiple spec at input.\n"
00348                   "Do not mix surface input types together\n");
00349       exit(1);
00350    }
00351 
00352       if (Spec->N_Surfs != 2) {
00353       SUMA_S_Err("2 surfaces expected.");
00354       exit(1);
00355    }
00356    
00357    SO1 = SUMA_Load_Spec_Surf(Spec, 0, ps->sv[0], 0);
00358    if (!SO1) {
00359          fprintf (SUMA_STDERR,"Error %s:\n"
00360                               "Failed to find surface\n"
00361                               "in spec file. \n",
00362                               FuncName );
00363          exit(1);
00364    }
00365    if (!SUMA_SurfaceMetrics(SO1, "EdgeList|MemberFace", NULL)) { SUMA_SL_Err("Failed to create edge list for SO1"); exit(1);  }
00366    if (Opt->fix_winding) {
00367             int orient, trouble;
00368             if (LocalHead) fprintf(SUMA_STDERR,"%s: Making sure S1 is consistently orientated\n", FuncName);
00369             if (!SUMA_MakeConsistent (SO1->FaceSetList, SO1->N_FaceSet, SO1->EL, Opt->debug, &trouble)) {
00370                SUMA_SL_Err("Failed in SUMA_MakeConsistent");
00371             }
00372             if (trouble && LocalHead) {
00373                fprintf(SUMA_STDERR,"%s: trouble value of %d from SUMA_MakeConsistent.\n"
00374                                     "Inconsistencies were found and corrected unless \n"
00375                                     "stderr output messages from SUMA_MakeConsistent\n"
00376                                     "indicate otherwise.\n", FuncName, trouble);
00377             }
00378             if (LocalHead) fprintf(SUMA_STDERR,"%s: Checking orientation.\n", FuncName);
00379             orient = SUMA_OrientTriangles (SO1->NodeList, SO1->N_Node, SO1->FaceSetList, SO1->N_FaceSet, 1, 0);
00380             if (orient < 0) { 
00381                /* flipping was done, dump the edge list since it is not automatically updated (should do that in function, just like in SUMA_MakeConsistent,  shame on you) */ 
00382                if (SO1->EL) SUMA_free_Edge_List(SO1->EL); SO1->EL = NULL; 
00383                if (!SUMA_SurfaceMetrics(SO1, "EdgeList", NULL)) { SUMA_SL_Err("Failed to create edge list for SO1"); exit(1);  }
00384                /* free normals, new ones needed (Normals should be flipped inside  of SUMA_OrientTriangles! (just like in SUMA_MakeConsistent) ) */
00385                if (SO1->NodeNormList) SUMA_free(SO1->NodeNormList); SO1->NodeNormList = NULL;
00386                if (SO1->FaceNormList) SUMA_free(SO1->FaceNormList); SO1->FaceNormList = NULL;
00387             }
00388             if (!orient) { fprintf(SUMA_STDERR,"Error %s:\nFailed in SUMA_OrientTriangles\n", FuncName); }
00389             if (LocalHead) {
00390                if (orient < 0) { SUMA_SL_Note("S1 was reoriented"); }
00391                else { SUMA_SL_Note("S1 was properly oriented"); }
00392             }
00393    }
00394   
00395    
00396    if (!SO1->NodeNormList || !SO1->FaceNormList) { SUMA_LH("Node Normals"); SUMA_RECOMPUTE_NORMALS(SO1); }
00397    if (Opt->NodeDbg >= SO1->N_Node) {
00398       SUMA_SL_Warn("node_debug index is larger than number of nodes in surface, ignoring -node_debug.");
00399       Opt->NodeDbg = -1;
00400    }
00401       
00402    SO2 = SUMA_Load_Spec_Surf(Spec, 1, ps->sv[1], 0);
00403    if (!SO2) {
00404          fprintf (SUMA_STDERR,"Error %s:\n"
00405                               "Failed to find surface\n"
00406                               "in spec file. \n",
00407                               FuncName );
00408          exit(1);
00409    }  
00410    if (!SUMA_SurfaceMetrics(SO2, "EdgeList|MemberFace", NULL)) { SUMA_SL_Err("Failed to create edge list for SO2"); exit(1);  }
00411    if (!SO2->NodeNormList || !SO2->FaceNormList) { SUMA_LH("Node Normals"); SUMA_RECOMPUTE_NORMALS(SO2); }
00412    
00413    if (LocalHead) { 
00414       SUMA_LH("Surf1");
00415       SUMA_Print_Surface_Object(SO1, NULL);
00416       SUMA_LH("Surf2");
00417       SUMA_Print_Surface_Object(SO2, NULL);
00418    }
00419    
00420    /* a select list of nodes? */
00421    nodeind = NULL; N_nodeind = 0;
00422    if (Opt->in_nodeindices) {
00423       im = mri_read_1D(Opt->in_nodeindices);
00424       if (!im) { SUMA_SL_Err("Failed to read 1D file of node indices"); exit(1);}
00425       far = MRI_FLOAT_PTR(im);
00426       N_nodeind = nvec = im->nx;
00427       ncol = im->ny;
00428       if (ncol != 1) { SUMA_SL_Err("More than one column in node index input file."); exit(1);}
00429       nodeind = (int *)SUMA_calloc(nvec, sizeof(int));
00430       if (!nodeind) { SUMA_SL_Crit("Failed to allocate"); exit(1); }
00431       for (i=0;i<nvec;++i) { 
00432          nodeind[i] = (int)far[i]; 
00433          if (nodeind[i] < 0 || nodeind[i] >= SO1->N_Node) {
00434             fprintf(SUMA_STDERR, "Error %s: A node index of %d was found in input file %s, entry %d.\n"
00435                                  "Acceptable indices are positive and less than %d\n", 
00436                                     FuncName, nodeind[i], Opt->in_nodeindices, i, SO1->N_Node);
00437             exit(1);
00438          }
00439       } 
00440       mri_free(im); im = NULL;   /* done with that baby */
00441    }
00442    
00443    /* a preset directions vector ?*/
00444    projdir = NULL; 
00445    if (Opt->in_1D) {
00446       im = mri_read_1D(Opt->in_1D);
00447       if (!im) { SUMA_SL_Err("Failed to read 1D file of projection directions"); exit(1);}
00448       far = MRI_FLOAT_PTR(im);
00449       if (im->ny != 3) { SUMA_SL_Err("Need three columns in projection directions file."); exit(1); }
00450       if (im->nx != SO1->N_Node) {
00451          fprintf(SUMA_STDERR, "Error %s: You must have a direction for each node in SO1.\n"
00452                               "%d directions found but SO1 has %d nodes.\n", FuncName, im->nx, SO1->N_Node);
00453          exit(1);
00454       }
00455 
00456       /* change to row major major and make it match nodeind */
00457       projdir = (float *)SUMA_calloc(SO1->N_Node*3, sizeof(float));
00458       if (!projdir) { SUMA_SL_Crit("Failed to allocate"); exit(1); }
00459       for (i=0; i<SO1->N_Node; ++i) {
00460          projdir[3*i  ] = far[i              ];
00461          projdir[3*i+1] = far[i+  SO1->N_Node];
00462          projdir[3*i+2] = far[i+2*SO1->N_Node];
00463       }
00464       mri_free(im); im = NULL;   /* done with that baby */
00465 
00466    }
00467    
00468    if (SO_name) {
00469       /* user is interpolating surface coords, check on other input insanity */
00470       if (nodeind) {
00471          fprintf(SUMA_STDERR, "Error %s: You cannot combine option -o_TYPE with -node_indices", FuncName);
00472          exit(1);
00473       }
00474       if (Opt->in_name) {
00475          fprintf(SUMA_STDERR, "Error %s: You cannot combine option -o_TYPE with -data", FuncName);
00476          exit(1);
00477       }
00478    } 
00479    /* a file containing data? */
00480    if (Opt->in_name) {
00481       im_data = mri_read_1D(Opt->in_name);
00482       if (!im_data) { SUMA_SL_Err("Failed to read 1D file of data"); exit(1);}
00483       far_data = MRI_FLOAT_PTR(im_data);
00484       nvec_data = im_data->nx;
00485       ncol_data = im_data->ny;
00486       if (nvec_data != SO2->N_Node) {
00487          SUMA_SL_Err("Your data file must have one row for each node in surface 2.\n"); exit(1);
00488       }
00489       d_order = SUMA_COLUMN_MAJOR;
00490    } else { 
00491       im_data = NULL;
00492       far_data = SO2->NodeList;
00493       nvec_data = SO2->N_Node;
00494       ncol_data = 3;
00495       d_order = SUMA_ROW_MAJOR;
00496    }
00497    
00498      
00499    SUMA_LH("Going for the mapping of SO1 --> SO2");
00500    M2M = SUMA_GetM2M_NN( SO1, SO2, nodeind, N_nodeind, projdir, 0, Opt->NodeDbg);
00501    
00502    /* Now show the mapping results for a debug node ? */
00503    if (Opt->NodeDbg >= 0) {
00504       char *s = NULL;
00505       s = SUMA_M2M_node_Info(M2M, Opt->NodeDbg);
00506       fprintf(SUMA_STDERR,"%s: Debug for node %d ([%f, %f, %f])of SO1:\n%s\n\n", 
00507                            FuncName, Opt->NodeDbg, 
00508                            SO1->NodeList[3*Opt->NodeDbg], SO1->NodeList[3*Opt->NodeDbg+1], SO1->NodeList[3*Opt->NodeDbg+2],
00509                            s); 
00510       SUMA_free(s); s = NULL;
00511    }
00512    
00513    /* Now please do the interpolation */
00514    if (Opt->Data) {
00515       if (Opt->NearestNode > 1) dt = SUMA_M2M_interpolate(M2M, far_data, ncol_data, nvec_data, d_order, 0 );
00516       else if (Opt->NearestNode == 1) dt = SUMA_M2M_interpolate(M2M, far_data, ncol_data, nvec_data, d_order, 1 );
00517       if (!dt) {
00518          SUMA_SL_Err("Failed to interpolate");
00519          exit(1);
00520       }
00521    }
00522    
00523    SUMA_LH("Forming the output");
00524    outptr = fopen(outname,"w");
00525    if (!outptr) {
00526       SUMA_SL_Err("Failed to open file for output.\n");
00527       exit(1);
00528    }
00529    
00530    /* first create the header of the output */
00531    SS = SUMA_StringAppend(NULL, NULL);
00532    SS = SUMA_StringAppend_va(SS, "#Mapping from nodes on surf 1 (S1) to nodes on surf 2 (S2)\n"
00533                                  "#  Surf 1 is labeled %s, idcode:%s\n"
00534                                  "#  Surf 2 is labeled %s, idcode:%s\n",
00535                                  SO1->Label, SO1->idcode_str, SO2->Label, SO2->idcode_str);
00536    icol = 0;
00537    SS = SUMA_StringAppend_va(SS, "#Col. %d:\n"
00538                                  "#     S1n (or nj): Index of node on S1\n"
00539                                  , icol); ++icol;
00540    if (Opt->NearestNode > 1) {
00541       SS = SUMA_StringAppend_va(SS, 
00542                                  "#Col. %d..%d:\n"
00543                                  "#     S2ne_S1n: Indices of %d nodes on S2 \n"
00544                                  "#     that are closest neighbors of nj.\n"
00545                                  "#     The first index is that of the node on S2 that is closest to nj.\n"
00546                                  "#     If -1 then thes values should be ignored. This happens when nj's projection failed.\n" 
00547                                  , icol, icol+Opt->NearestNode-1, Opt->NearestNode); icol += Opt->NearestNode;
00548       SS = SUMA_StringAppend_va(SS, 
00549                                  "#Col. %d..%d:\n"
00550                                  "#     S2we_S1n: Weights assigned to nodes on surf 2 (S2) \n"
00551                                  "#     that are closest neighbors of nj.\n"
00552                                  , icol, icol+Opt->NearestNode-1, Opt->NearestNode); icol += Opt->NearestNode;
00553    } else if (Opt->NearestNode == 1) {
00554       SS = SUMA_StringAppend_va(SS, 
00555                                  "#Col. %d:\n"
00556                                  "#     S2ne_S1n: Index of the node on S2 (label:%s idcode:%s)\n"
00557                                  "#     that is the closest neighbor of nj.\n"
00558                                  "#     If -1 then this value should be ignored. This happens when nj's projection failed.\n" 
00559                                  , icol, SO2->Label, SO2->idcode_str); ++icol;
00560    }
00561    if (Opt->NearestTriangle) { 
00562       SS = SUMA_StringAppend_va(SS, 
00563                                  "#Col. %d:\n"
00564                                  "#     S2t_S1n: Index of the triangle on S2 that hosts node nj on S1.\n"
00565                                  "#     In other words, nj's closest projection onto S2 falls on \n"
00566                                  "#     triangle S2t_S1n\n"
00567                                  "#     If -1 then this value should be ignored. This happens when nj's projection failed.\n" 
00568                                  , icol); ++icol; 
00569    }
00570    if (Opt->ProjectionOnMesh) { 
00571       SS = SUMA_StringAppend_va(SS, 
00572                                  "#Col. %d..%d:\n"
00573                                  "#     S2p_S1n: Coordinates of projection of nj onto S2\n"
00574                                  , icol, icol+2); icol += 3; 
00575    }
00576    if (Opt->DistanceToMesh) {
00577       SS = SUMA_StringAppend_va(SS, 
00578                                  "#Col. %d:\n"
00579                                  "#     Closest distance from nj to S2\n"
00580                                  , icol); ++icol;
00581    }
00582    if (Opt->Data) {
00583       if (!Opt->in_name) {
00584          SS = SUMA_StringAppend_va(SS, 
00585                                     "#Col. %d..%d:\n"
00586                                     "#     Interpolation using XYZ coordinates of nodes on S2 that neighbor nj\n"
00587                                     "#     (same as coordinates of node's projection onto triangle in S2, if using \n"
00588                                     "       barycentric interpolation)\n"
00589                                     , icol, icol+2); icol += 3; 
00590       } else {
00591          SS = SUMA_StringAppend_va(SS, 
00592                                     "#Col. %d..%d:\n"
00593                                     "#     Interpolation of data at nodes on S2 that neighbor nj\n"
00594                                     "#     Data obtained from %s\n"
00595                                     , icol, icol+ncol_data-1, Opt->in_name);  icol += ncol_data;
00596       }
00597    } 
00598    s = SUMA_HistString("SurfToSurf", argc, argv, NULL);
00599    SS = SUMA_StringAppend_va(SS, 
00600                                 "#History:\n"
00601                                 "#%s\n", s); SUMA_free(s); s = NULL;
00602    SUMA_SS2S(SS,s);
00603    fprintf(outptr,"%s\n",s); SUMA_free(s); s = NULL;
00604    
00605    /* put headers atop columns */
00606    Nchar = 6; /* if you change this number you'll need to fix the formats below */
00607    for (i=0; i<icol; ++i) { sprintf(sbuf,"#%s", MV_format_fval2(i, Nchar -1)); fprintf(outptr,"%6s   ", sbuf); }
00608    fprintf(outptr,"\n");
00609    
00610    /* Now put in the values, make sure you parallel columns above! */
00611    for (i=0; i<M2M->M1Nn; ++i) {
00612       fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M1n[i], Nchar));
00613       if (Opt->NearestNode > 0) {
00614          for (j=0; j<Opt->NearestNode; ++j) { 
00615             if (j < M2M->M2Nne_M1n[i]) fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2ne_M1n[i][j], Nchar)); 
00616             else fprintf(outptr,"%6s   ", "-1"); 
00617          } /* Neighboring nodes */
00618       } 
00619       if (Opt->NearestNode > 1) { /* add the weights */
00620          for (j=0; j<Opt->NearestNode; ++j) { 
00621             if (j < M2M->M2Nne_M1n[i]) fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2we_M1n[i][j], Nchar)); 
00622             else fprintf(outptr,"%6s   ", "0.0"); 
00623          } 
00624       }
00625       if (Opt->NearestTriangle) {
00626          fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2t_M1n[i], Nchar)); 
00627       }
00628       if (Opt->ProjectionOnMesh) {
00629          fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2p_M1n[3*i], Nchar));
00630          fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2p_M1n[3*i+1], Nchar));
00631          fprintf(outptr,"%6s   ", MV_format_fval2(M2M->M2p_M1n[3*i+2], Nchar)); 
00632       }
00633       if (Opt->DistanceToMesh) { 
00634          fprintf(outptr,"%6s   ", MV_format_fval2(M2M->PD[i], Nchar)); 
00635       }
00636       if (Opt->Data) {
00637          if (!Opt->in_name) {
00638             fprintf(outptr,"%6s   ", MV_format_fval2(dt[3*i], Nchar));
00639             fprintf(outptr,"%6s   ", MV_format_fval2(dt[3*i+1], Nchar));
00640             fprintf(outptr,"%6s   ", MV_format_fval2(dt[3*i+2], Nchar));
00641          } else { /* Column major business */
00642             for (j=0; j<ncol_data; ++j) { fprintf(outptr,"%6s   ", MV_format_fval2(dt[i+j*M2M->M1Nn], Nchar)); }
00643          }
00644       }
00645       fprintf(outptr,"\n");
00646    }
00647    
00648    /* do they want an output surface ? */
00649    if (SO_name) {
00650       float *tmpfv = NULL;
00651       SUMA_LH("Writing surface");
00652       tmpfv = SO1->NodeList;
00653       SO1->NodeList = dt;
00654       if (!SUMA_Save_Surface_Object (SO_name, SO1, ps->o_FT[0], ps->o_FF[0], NULL)) {
00655          fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
00656          exit (1);
00657       }
00658       SO1->NodeList = tmpfv; tmpfv = NULL;
00659    }
00660    
00661    if (projdir) SUMA_free(projdir); projdir = NULL;
00662    if (SO_name) SUMA_free(SO_name); SO_name = NULL;   
00663    if (outptr) fclose(outptr); outptr = NULL;
00664    if (dt) SUMA_free(dt); dt = NULL;
00665    if (s) SUMA_free(s); s = NULL;
00666    if (im_data) mri_free(im_data); im_data = NULL;   /* done with the data */
00667    if (nodeind) SUMA_free(nodeind); nodeind = NULL;
00668    if (M2M) M2M = SUMA_FreeM2M(M2M);
00669    if (SO1) SUMA_Free_Surface_Object(SO1); SO1 = NULL;
00670    if (SO2) SUMA_Free_Surface_Object(SO2); SO2 = NULL;
00671    if (Spec) SUMA_free(Spec); Spec = NULL;
00672    if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
00673    if (Opt) Opt = SUMA_Free_Generic_Prog_Options_Struct(Opt);
00674    if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
00675    exit(0);
00676    
00677 } 
 

Powered by Plone

This site conforms to the following standards: