1 #include "SUMA_suma.h"
2
usage_SurfToSurf(SUMA_GENERIC_ARGV_PARSE * ps,int detail)3 void usage_SurfToSurf (SUMA_GENERIC_ARGV_PARSE *ps, int detail)
4 {
5 static char FuncName[]={"usage_SurfToSurf"};
6 char * s = NULL, *sio=NULL, *st = NULL, *sts = NULL;
7 int i;
8 s = SUMA_help_basics();
9 sio = SUMA_help_IO_Args(ps);
10 printf (
11 "\n"
12 "Usage: SurfToSurf <-i_TYPE S1> [<-sv SV1>]\n"
13 " <-i_TYPE S2> [<-sv SV1>]\n"
14 " [<-prefix PREFIX>]\n"
15 " [<-output_params PARAM_LIST>]\n"
16 " [<-node_indices NODE_INDICES>]\n"
17 " [<-proj_dir PROJ_DIR>]\n"
18 " [<-data DATA>]\n"
19 " [<-node_debug NODE>]\n"
20 " [<-debug DBG_LEVEL>]\n"
21 " [-make_consistent]\n"
22 " [<-dset DSET>]\n"
23 " [<-mapfile MAP_INFO>]\n"
24 " \n"
25 " This program is used to interpolate data from one surface (S2)\n"
26 " to another (S1), assuming the surfaces are quite similar in\n"
27 " shape but having different meshes (non-isotopic).\n"
28 " This is done by projecting each node (nj) of S1 along the normal\n"
29 " at nj and finding the closest triangle t of S2 that is intersected\n"
30 " by this projection. Projection is actually bidirectional.\n"
31 " If such a triangle t is found, the nodes (of S2) forming it are \n"
32 " considered to be the neighbors of nj.\n"
33 " Values (arbitrary data, or coordinates) at these neighboring nodes\n"
34 " are then transfered to nj using barycentric interpolation or \n"
35 " nearest-node interpolation.\n"
36 " Nodes whose projections fail to intersect triangles in S2 are given\n"
37 " nonsensical values of -1 and 0.0 in the output.\n"
38 "\n"
39 " Mandatory input:\n"
40 " Two surfaces are required at input. See -i_TYPE options\n"
41 " below for more information. \n"
42 "\n"
43 " Optional input:\n"
44 " -prefix PREFIX: Specify the prefix of the output file.\n"
45 " The output file is in 1D format at the moment.\n"
46 " Default is SurfToSurf\n"
47 " -output_params PARAM_LIST: Specify the list of mapping\n"
48 " parameters to include in output\n"
49 " PARAM_LIST can have any or all of the following:\n"
50 " NearestTriangleNodes: Use Barycentric interpolation (default)\n"
51 " and output indices of 3 nodes from S2\n"
52 " that neighbor nj of S1\n"
53 " NearestNode: Use only the closest node from S2 (of the three \n"
54 " closest neighbors) to nj of S1 for interpolation\n"
55 " and output the index of that closest node.\n"
56 " NearestTriangle: Output index of triangle t from S2 that\n"
57 " is the closest to nj along its projection\n"
58 " direction. \n"
59 " DistanceToSurf: Output distance (signed) from nj, along \n"
60 " projection direction to S2.\n"
61 " This is the parameter output by the precursor\n"
62 " program CompareSurfaces\n"
63 " ProjectionOnSurf: Output coordinates of projection of nj onto \n"
64 " triangle t of S2.\n"
65 " NearestNodeCoords: X Y Z coordinates of closest node on S2\n"
66 " Data: Output the data from S2, interpolated onto S1\n"
67 " If no data is specified via the -data option, then\n"
68 " the XYZ coordinates of SO2's nodes are considered\n"
69 " the data.\n"
70 " -data DATA: 1D file containing data to be interpolated.\n"
71 " Each row i contains data for node i of S2.\n"
72 " You must have one row for each node making up S2.\n"
73 " In other terms, if S2 has N nodes, you need N rows\n"
74 " in DATA. \n"
75 " Each column of DATA is processed separately (think\n"
76 " sub-bricks, and spatial interpolation).\n"
77 " You can use [] selectors to choose a subset \n"
78 " of columns.\n"
79 " If -data option is not specified and Data is in PARAM_LIST\n"
80 " then the XYZ coordinates of SO2's nodes are the data.\n"
81 " -dset DSET: Treat like -data, but works best with datasets, preserving\n"
82 " header information in the output.\n"
83 " -dset and -data are mutually exclusive.\n"
84 " Also, -dset and parameter Data cannot be mixed.\n"
85 " -node_indices NODE_INDICES: 1D file containing the indices of S1\n"
86 " to consider. The default is all of the\n"
87 " nodes in S1. Only one column of values is\n"
88 " allowed here, use [] selectors to choose\n"
89 " the column of node indices if NODE_INDICES\n"
90 " has multiple columns in it.\n"
91 " -proj_dir PROJ_DIR: 1D file containing projection directions to use\n"
92 " instead of the node normals of S1.\n"
93 " Each row should contain one direction for each\n"
94 " of the nodes forming S1.\n"
95 " -closest_possible OO: Flag allowing the substitution of the projection\n"
96 " result with the closest node that could be found\n"
97 " along any direction.\n"
98 " 0: Don't do that, direction results only.\n"
99 " 1: Use closest node if projection fails to hit target\n"
100 " 2: Use closest node if it is at a closer distance.\n"
101 " 3: Use closest and don't bother with projections.\n"
102 " -make_consistent: Force a consistency check and correct triangle \n"
103 " orientation of S1 if needed. Triangles are also\n"
104 " oriented such that the majority of normals point\n"
105 " away from center of surface.\n"
106 " The program might not succeed in repairing some\n"
107 " meshes with inconsistent orientation.\n"
108 " -mapfile MAP_INFO: Use the mapping from S2 to S1 that is stored in\n"
109 " MAP_INFO. MAP_INFO is a file containing the mapping\n"
110 " parameters between surfaces S2 and S1. \n"
111 " It is generated automatically by SurfToSurf when \n"
112 " -mapfile is not used, and saved under PREFIX.niml.M2M.\n"
113 " Reusing the MAP_INFO file allows for faster execution\n"
114 " of SurfToSurf the next time around, assuming of course\n"
115 " that the two surfaces involved are the same, and that \n"
116 " only the input data differs.\n"
117 " MAP_INFO is also generated by MapIcosahedron.\n"
118 "\n"
119 "%s"
120 "%s"
121 "\n", sio, s);
122 SUMA_free(s); s = NULL; SUMA_free(st); st = NULL; SUMA_free(sio); sio = NULL;
123 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
124 printf(" Ziad S. Saad SSCC/NIMH/NIH saadz@mail.nih.gov \n"
125 " Shruti Japee LBC/NIMH/NIH shruti@codon.nih.gov \n");
126 exit(0);
127 }
128
SUMA_SurfToSurf_ParseInput(char * argv[],int argc,SUMA_GENERIC_ARGV_PARSE * ps)129 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_SurfToSurf_ParseInput(
130 char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
131 {
132 static char FuncName[]={"SUMA_BrainWrap_ParseInput"};
133 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
134 int kar;
135 SUMA_Boolean brk, accepting_out;
136
137 SUMA_Boolean LocalHead = NOPE;
138
139 SUMA_ENTRY;
140
141 Opt = SUMA_Alloc_Generic_Prog_Options_Struct();
142 kar = 1;
143 brk = NOPE;
144 Opt->in_1D = NULL;
145 Opt->NodeDbg = -1;
146 Opt->debug = 0;
147 Opt->NearestNode = 0;
148 Opt->NearestTriangle = 0;
149 Opt->DistanceToMesh = 0;
150 Opt->ProjectionOnMesh = 0;
151 Opt->NearestNodeCoords = 0;
152 Opt->Data = 0;
153 Opt->in_name = NULL;
154 Opt->out_prefix = NULL;
155 Opt->fix_winding = 0;
156 Opt->iopt = 0;
157 Opt->oform = SUMA_NO_DSET_FORMAT;
158 accepting_out = NOPE;
159 while (kar < argc) { /* loop accross command ine options */
160 /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
161
162 if (!brk && accepting_out) {
163 /* make sure you have not begun with new options */
164 if (*(argv[kar]) == '-') accepting_out = NOPE;
165 }
166
167 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
168 usage_SurfToSurf(ps, strlen(argv[kar]) > 3 ? 2:1);
169 exit (0);
170 }
171
172 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
173
174 if (!brk && (strcmp(argv[kar], "-debug") == 0))
175 {
176 if (kar+1 >= argc)
177 {
178 fprintf (SUMA_STDERR, "need a number after -debug \n");
179 exit (1);
180 }
181
182 Opt->debug = atoi(argv[++kar]);
183 brk = YUP;
184 }
185
186 if (!brk && (strcmp(argv[kar], "-node_debug") == 0))
187 {
188 if (kar+1 >= argc)
189 {
190 fprintf (SUMA_STDERR, "need a number after -node_debug \n");
191 exit (1);
192 }
193
194 Opt->NodeDbg = atoi(argv[++kar]);
195 brk = YUP;
196 }
197
198 if (!brk && (strcmp(argv[kar], "-node_indices") == 0))
199 {
200 if (kar+1 >= argc)
201 {
202 fprintf (SUMA_STDERR, "need a parameter after -node_indices \n");
203 exit (1);
204 }
205
206 Opt->in_nodeindices = argv[++kar];
207 brk = YUP;
208 }
209
210 if (!brk && (strcmp(argv[kar], "-closest_possible") == 0))
211 {
212 if (kar+1 >= argc)
213 {
214 fprintf (SUMA_STDERR, "need a number after -closest_possible \n");
215 exit (1);
216 }
217
218 Opt->iopt = atoi(argv[++kar]);
219 if (Opt->iopt != 0 && Opt->iopt != 1 && Opt->iopt != 2 &&
220 Opt->iopt != 3) {
221 SUMA_S_Errv("Must choose from 0, 1, 2, or 3 for -closest_possible."
222 " Have %d\n",
223 Opt->iopt);
224 exit (1);
225 }
226 brk = YUP;
227 }
228
229 if (!brk && (strcmp(argv[kar], "-output_params") == 0))
230 {
231 if (kar+1 >= argc)
232 {
233 fprintf (SUMA_STDERR,
234 "need at least one parameter after output_params \n");
235 exit (1);
236 }
237
238 accepting_out = YUP;
239 brk = YUP;
240 }
241
242 if (!brk && accepting_out && (strcmp(argv[kar], "NearestNode") == 0)) {
243 if (Opt->NearestNode < 1) Opt->NearestNode = 1;
244 brk = YUP;
245 }
246
247 if (!brk && accepting_out &&
248 (strcmp(argv[kar], "NearestTriangleNodes") == 0)) {
249 if (Opt->NearestNode < 3) Opt->NearestNode = 3;
250 brk = YUP;
251 }
252
253 if (!brk && accepting_out && (strcmp(argv[kar], "NearestTriangle") == 0)) {
254 if (Opt->NearestTriangle < 1) Opt->NearestTriangle = 1;
255 brk = YUP;
256 }
257
258 if (!brk && accepting_out && (strcmp(argv[kar], "DistanceToSurf") == 0)) {
259 Opt->DistanceToMesh = 1;
260 brk = YUP;
261 }
262
263 if (!brk && accepting_out &&
264 (strcmp(argv[kar], "ProjectionOnSurf") == 0)) {
265 Opt->ProjectionOnMesh = 1;
266 brk = YUP;
267 }
268
269 if (!brk && accepting_out &&
270 (strcmp(argv[kar], "NearestNodeCoords") == 0)) {
271 Opt->NearestNodeCoords = 1;
272 brk = YUP;
273 }
274
275 if (!brk && accepting_out &&
276 (strcmp(argv[kar], "Data") == 0)) {
277 if (Opt->Data < 0) {
278 fprintf (SUMA_STDERR,
279 "Cannot mix parameter Data with -dset option \n");
280 exit (1);
281 }
282 Opt->Data = 1;
283 brk = YUP;
284 }
285
286 if (!brk && (strcmp(argv[kar], "-data") == 0))
287 {
288 if (kar+1 >= argc)
289 {
290 fprintf (SUMA_STDERR, "need a name after -data \n");
291 exit (1);
292 }
293 ++kar;
294 if (strcmp(argv[kar],"_XYZ_") == 0) {
295 /* default Opt->in_name = NULL*/
296 if (Opt->in_name) {
297 SUMA_SL_Err("Input already specified."
298 "Do not mix -data and -dset");
299 exit (1);
300 }
301 } else {
302 Opt->in_name = SUMA_copy_string(argv[kar]);
303 }
304 Opt->Data = 1;
305 brk = YUP;
306 }
307
308 if (!brk && (strcmp(argv[kar], "-dset") == 0))
309 {
310 if (kar+1 >= argc)
311 {
312 fprintf (SUMA_STDERR, "need a name after -dset \n");
313 exit (1);
314 }
315 ++kar;
316 if (strcmp(argv[kar],"_XYZ_") == 0 || Opt->Data > 0) {
317 /* default Opt->in_name = NULL*/
318 if (Opt->in_name || Opt->Data > 0) {
319 SUMA_SL_Err("Input already specified."
320 "Do not mix -data and -dset."
321 "Or use parameter DATA with -dset");
322 exit (1);
323 }
324 } else {
325 Opt->in_name = SUMA_copy_string(argv[kar]);
326 }
327 Opt->Data = -1;
328 brk = YUP;
329 }
330
331 if (!brk && (strcmp(argv[kar], "-prefix") == 0))
332 {
333 if (kar+1 >= argc)
334 {
335 fprintf (SUMA_STDERR, "need a name after -prefix \n");
336 exit (1);
337 }
338
339 Opt->out_prefix = SUMA_RemoveDsetExtension_eng(argv[++kar],
340 &(Opt->oform));
341 brk = YUP;
342 }
343
344 if (!brk && (strcmp(argv[kar], "-mapfile") == 0))
345 {
346 if (kar+1 >= argc)
347 {
348 fprintf (SUMA_STDERR, "need a name after -mapfile \n");
349 exit (1);
350 }
351
352 Opt->s = SUMA_Extension(argv[++kar],".niml.M2M", NOPE);
353 if (!SUMA_filexists(Opt->s)) {
354 SUMA_S_Errv("File %s not found\n"
355 , Opt->s);
356 exit(1);
357 }
358 brk = YUP;
359 }
360
361 if (!brk && (strcmp(argv[kar], "-proj_dir") == 0))
362 {
363 if (kar+1 >= argc)
364 {
365 fprintf (SUMA_STDERR, "need a name after -proj_dir \n");
366 exit (1);
367 }
368
369 Opt->in_1D = argv[++kar];
370 brk = YUP;
371 }
372 if (!brk && (strcmp(argv[kar], "-make_consistent") == 0))
373 {
374 Opt->fix_winding = 1;
375 brk = YUP;
376 }
377 if (!brk && !ps->arg_checked[kar]) {
378 fprintf (SUMA_STDERR,
379 "Error %s:\n"
380 "Option %s not understood. Try -help for usage\n",
381 FuncName, argv[kar]);
382 suggest_best_prog_option(argv[0], argv[kar]);
383 exit (1);
384 } else {
385 brk = NOPE;
386 kar ++;
387 }
388 }
389
390 /* set default for NearestNode if nothing has been set */
391 if (Opt->NearestNode < 1) Opt->NearestNode = 3;
392
393 if (!Opt->out_prefix) Opt->out_prefix = SUMA_copy_string("SurfToSurf");
394 if (Opt->in_1D && Opt->s) {
395 SUMA_S_Err("Cannot use -proj_dir along with -mapfile");
396 exit(1);
397 }
398 SUMA_RETURN(Opt);
399 }
400
main(int argc,char * argv[])401 int main (int argc,char *argv[])
402 {/* Main */
403 static char FuncName[]={"SurfToSurf"};
404 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
405 SUMA_GENERIC_ARGV_PARSE *ps=NULL;
406 SUMA_SurfaceObject *SO1=NULL, *SO2 = NULL;
407 SUMA_SurfSpecFile *Spec = NULL;
408 SUMA_M2M_STRUCT *M2M = NULL;
409 int N_Spec=0, *nodeind = NULL, N_nodeind, icol, i, j;
410 MRI_IMAGE *im = NULL, *im_data=NULL;
411 int nvec=0, ncol=0, nvec_data=0, ncol_data=0, Nchar=0;
412 float *far = NULL, *far_data=NULL, *dt = NULL, *projdir=NULL;
413 char *outname = NULL, *s=NULL, sbuf[100];
414 void *SO_name = NULL;
415 FILE *outptr=NULL;
416 SUMA_Boolean exists = NOPE;
417 SUMA_INDEXING_ORDER d_order = SUMA_NO_ORDER;
418 SUMA_STRING *SS=NULL;
419 SUMA_Boolean LocalHead = NOPE;
420
421 SUMA_STANDALONE_INIT;
422 SUMA_mainENTRY;
423
424 /* Allocate space for DO structure */
425 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
426 ps = SUMA_Parse_IO_Args(argc, argv, "-i;-t;-spec;-s;-sv;-o;");
427
428 Opt = SUMA_SurfToSurf_ParseInput (argv, argc, ps);
429 if (argc < 2) {
430 SUMA_S_Err("Too few options");
431 usage_SurfToSurf(ps, 0);
432 exit (1);
433 }
434
435
436 /* if output surface requested, check on pre-existing file */
437 if (ps->o_N_surfnames) {
438 SO_name = SUMA_Prefix2SurfaceName(ps->o_surfnames[0],
439 NULL, NULL, ps->o_FT[0], &exists);
440 if (exists) {
441 fprintf(SUMA_STDERR,
442 "Error %s:\nOutput file(s) %s* on disk.\n"
443 "Will not overwrite.\n", FuncName, ps->o_surfnames[0]);
444 exit(1);
445 }
446 }
447
448 if (Opt->debug > 2) LocalHead = YUP;
449 outname = SUMA_append_extension(Opt->out_prefix,".1D");
450 if (SUMA_filexists(outname) && !THD_ok_overwrite()) {
451 fprintf(SUMA_STDERR,"Output file %s exists.\n", outname);
452 exit(1);
453 }
454
455 /* Load the surfaces from command line*/
456 Spec = SUMA_IO_args_2_spec(ps, &N_Spec);
457 if (N_Spec != 1) {
458 SUMA_S_Err( "Multiple spec at input.\n"
459 "Do not mix surface input types together\n");
460 exit(1);
461 }
462
463 if (Spec->N_Surfs != 2) {
464 SUMA_S_Err("2 surfaces expected.");
465 exit(1);
466 }
467
468 SO1 = SUMA_Load_Spec_Surf(Spec, 0, ps->sv[0], 0);
469 if (!SO1) {
470 fprintf (SUMA_STDERR,"Error %s:\n"
471 "Failed to find surface\n"
472 "in spec file. \n",
473 FuncName );
474 exit(1);
475 }
476 if (!SUMA_SurfaceMetrics(SO1, "EdgeList|MemberFace", NULL)) {
477 SUMA_SL_Err("Failed to create edge list for SO1");
478 exit(1);
479 }
480 if (Opt->fix_winding) {
481 int orient, trouble;
482 if (LocalHead)
483 fprintf(SUMA_STDERR,
484 "%s: Making sure S1 is consistently orientated\n", FuncName);
485 if (!SUMA_MakeConsistent (SO1->FaceSetList, SO1->N_FaceSet,
486 SO1->EL, Opt->debug, &trouble)) {
487 SUMA_SL_Err("Failed in SUMA_MakeConsistent");
488 }
489 if (trouble && LocalHead) {
490 fprintf(SUMA_STDERR,
491 "%s: trouble value of %d from SUMA_MakeConsistent.\n"
492 "Inconsistencies were found and corrected unless \n"
493 "stderr output messages from SUMA_MakeConsistent\n"
494 "indicate otherwise.\n", FuncName, trouble);
495 }
496 if (LocalHead)
497 fprintf(SUMA_STDERR,"%s: Checking orientation.\n", FuncName);
498 orient = SUMA_OrientTriangles (SO1->NodeList, SO1->N_Node,
499 SO1->FaceSetList, SO1->N_FaceSet,
500 1, 0, NULL, NULL);
501 if (orient < 0) {
502 /* flipping was done, dump the edge list since it is not
503 automatically updated (should do that in function,
504 just like in SUMA_MakeConsistent, shame on you)
505
506 If you revisit this section, use the newer:
507 SUMA_OrientSOTriangles
508 */
509 if (SO1->EL) SUMA_free_Edge_List(SO1->EL); SO1->EL = NULL;
510 if (!SUMA_SurfaceMetrics(SO1, "EdgeList", NULL)) {
511 SUMA_SL_Err("Failed to create edge list for SO1"); exit(1);
512 }
513 /* free normals, new ones needed (Normals should be flipped inside of
514 SUMA_OrientTriangles! (just like in SUMA_MakeConsistent) ) */
515 if (SO1->NodeNormList) SUMA_free(SO1->NodeNormList);
516 SO1->NodeNormList = NULL;
517 if (SO1->FaceNormList) SUMA_free(SO1->FaceNormList);
518 SO1->FaceNormList = NULL;
519 }
520 if (!orient) {
521 fprintf(SUMA_STDERR,
522 "Error %s:\nFailed in SUMA_OrientTriangles\n", FuncName);
523 }
524 if (LocalHead) {
525 if (orient < 0) { SUMA_SL_Note("S1 was reoriented"); }
526 else { SUMA_SL_Note("S1 was properly oriented"); }
527 }
528 }
529
530
531 if (!SO1->NodeNormList || !SO1->FaceNormList) {
532 SUMA_LH("Node Normals"); SUMA_RECOMPUTE_NORMALS(SO1);
533 }
534 if (Opt->NodeDbg >= SO1->N_Node) {
535 SUMA_SL_Warn( "node_debug index is larger than number "
536 "of nodes in surface, ignoring -node_debug.");
537 Opt->NodeDbg = -1;
538 }
539
540 SO2 = SUMA_Load_Spec_Surf(Spec, 1, ps->sv[1], 0);
541 if (!SO2) {
542 fprintf (SUMA_STDERR,"Error %s:\n"
543 "Failed to find surface\n"
544 "in spec file. \n",
545 FuncName );
546 exit(1);
547 }
548 if (!SUMA_SurfaceMetrics(SO2, "EdgeList|MemberFace", NULL)) {
549 SUMA_SL_Err("Failed to create edge list for SO2"); exit(1);
550 }
551 if (!SO2->NodeNormList || !SO2->FaceNormList) {
552 SUMA_LH("Node Normals"); SUMA_RECOMPUTE_NORMALS(SO2);
553 }
554
555 if (LocalHead) {
556 SUMA_LH("Surf1");
557 SUMA_Print_Surface_Object(SO1, NULL);
558 SUMA_LH("Surf2");
559 SUMA_Print_Surface_Object(SO2, NULL);
560 }
561
562 /* a select list of nodes? */
563 nodeind = NULL; N_nodeind = 0;
564 if (Opt->in_nodeindices) {
565 im = mri_read_1D(Opt->in_nodeindices);
566 if (!im) { SUMA_SL_Err("Failed to read 1D file of node indices"); exit(1);}
567 far = MRI_FLOAT_PTR(im);
568 N_nodeind = nvec = im->nx;
569 ncol = im->ny;
570 if (ncol != 1) {
571 SUMA_SL_Err("More than one column in node index input file."); exit(1);
572 }
573 nodeind = (int *)SUMA_calloc(nvec, sizeof(int));
574 if (!nodeind) { SUMA_SL_Crit("Failed to allocate"); exit(1); }
575 for (i=0;i<nvec;++i) {
576 nodeind[i] = (int)far[i];
577 if (nodeind[i] < 0 || nodeind[i] >= SO1->N_Node) {
578 fprintf(SUMA_STDERR,
579 "Error %s:\n"
580 "A node index of %d was found in input file %s, entry %d.\n"
581 "Acceptable indices are positive and less than %d\n",
582 FuncName, nodeind[i], Opt->in_nodeindices, i, SO1->N_Node);
583 exit(1);
584 }
585 }
586 mri_free(im); im = NULL; /* done with that baby */
587 }
588
589 /* a preset directions vector ?*/
590 projdir = NULL;
591 if (Opt->in_1D) {
592 im = mri_read_1D(Opt->in_1D);
593 if (!im) {
594 SUMA_SL_Err("Failed to read 1D file of projection directions"); exit(1);
595 }
596 far = MRI_FLOAT_PTR(im);
597 if (im->ny != 3) {
598 SUMA_SL_Err("Need three columns in projection directions file.");
599 exit(1);
600 }
601 if (im->nx != SO1->N_Node) {
602 fprintf(SUMA_STDERR,
603 "Error %s: You must have a direction for each node in SO1.\n"
604 "%d directions found but SO1 has %d nodes.\n",
605 FuncName, im->nx, SO1->N_Node);
606 exit(1);
607 }
608
609 /* change to row major major and make it match nodeind */
610 projdir = (float *)SUMA_calloc(SO1->N_Node*3, sizeof(float));
611 if (!projdir) { SUMA_SL_Crit("Failed to allocate"); exit(1); }
612 for (i=0; i<SO1->N_Node; ++i) {
613 projdir[3*i ] = far[i ];
614 projdir[3*i+1] = far[i+ SO1->N_Node];
615 projdir[3*i+2] = far[i+2*SO1->N_Node];
616 }
617 mri_free(im); im = NULL; /* done with that baby */
618
619 }
620
621 if (SO_name) {
622 /* user is interpolating surface coords, check on other input insanity */
623 if (nodeind) {
624 fprintf( SUMA_STDERR,
625 "Error %s: You cannot combine "
626 "option -o_TYPE with -node_indices", FuncName);
627 exit(1);
628 }
629 if (Opt->in_name) {
630 fprintf(SUMA_STDERR,
631 "Error %s: You cannot combine option -o_TYPE with -data",
632 FuncName);
633 exit(1);
634 }
635 }
636 /* a 1D file containing data, or Data parameter (for XYZ)? */
637 if (Opt->Data > 0) {
638 if (Opt->in_name) {
639 /* When you are ready to work with dsets, you should
640 checkout the function morphDsetToStd. It uses M2M */
641 im_data = mri_read_1D(Opt->in_name);
642 if (!im_data) {
643 SUMA_SL_Err("Failed to read 1D file of data"); exit(1);}
644 far_data = MRI_FLOAT_PTR(im_data);
645 nvec_data = im_data->nx;
646 ncol_data = im_data->ny;
647 if (nvec_data != SO2->N_Node) {
648 SUMA_SL_Err("Your data file must have one row "
649 "for each node in surface 2.\n"); exit(1);
650 }
651 d_order = SUMA_COLUMN_MAJOR;
652 } else {
653 im_data = NULL;
654 far_data = SO2->NodeList;
655 nvec_data = SO2->N_Node;
656 ncol_data = 3;
657 d_order = SUMA_ROW_MAJOR;
658 }
659 } else {
660 /* just -dset */
661 }
662
663
664
665
666 if (!Opt->s) {
667 SUMA_LH("Going for the mapping of SO1 --> SO2");
668 M2M = SUMA_GetM2M_NN( SO1, SO2, nodeind, N_nodeind,
669 projdir, 0, Opt->NodeDbg, Opt->iopt);
670 SUMA_S_Notev("Saving M2M into %s\n\n",
671 Opt->out_prefix);
672 if (!(SUMA_Save_M2M(Opt->out_prefix, M2M))) {
673 SUMA_S_Err("Failed to save M2M");
674 exit(1);
675 }
676 } else {
677 SUMA_S_Notev("Reusing mapping of SO1 --> SO2 from %s\n\n",
678 Opt->s);
679 if (!(M2M = SUMA_Load_M2M(Opt->s))) {
680 SUMA_S_Errv("Failed to load %s\n", Opt->s);
681 exit(1);
682 }
683 }
684
685 /* Now show the mapping results for a debug node ? */
686 if (Opt->NodeDbg >= 0) {
687 char *s = NULL;
688 s = SUMA_M2M_node_Info(M2M, Opt->NodeDbg);
689 fprintf(SUMA_STDERR,"%s: Debug for node %d ([%f, %f, %f])of SO1:\n%s\n\n",
690 FuncName, Opt->NodeDbg,
691 SO1->NodeList[3*Opt->NodeDbg],
692 SO1->NodeList[3*Opt->NodeDbg+1],
693 SO1->NodeList[3*Opt->NodeDbg+2],
694 s);
695 SUMA_free(s); s = NULL;
696 }
697
698 /* Now please do the interpolation */
699 if (Opt->Data > 0) {
700 if (Opt->NearestNode > 1)
701 dt = SUMA_M2M_interpolate( M2M, far_data, ncol_data,
702 nvec_data, d_order, 0 );
703 else if (Opt->NearestNode == 1)
704 dt = SUMA_M2M_interpolate( M2M, far_data, ncol_data,
705 nvec_data, d_order, 1 );
706 if (!dt) {
707 SUMA_SL_Err("Failed to interpolate");
708 exit(1);
709 }
710 } else if (Opt->Data < 0) {
711 SUMA_DSET *dset=NULL, *dseto=NULL;
712 char *oname=NULL, *uname=NULL, *s1=NULL, *s2=NULL;
713 int iform=SUMA_NO_DSET_FORMAT;
714 if (Opt->NodeDbg>= 0) {
715 SUMA_S_Notev("Processing dset %s\n", Opt->in_name);
716 }
717 iform = SUMA_NO_DSET_FORMAT;
718 if (!(dset = SUMA_LoadDset_s (Opt->in_name, &iform, 0))) {
719 SUMA_S_Errv("Failed to load %s\n", Opt->in_name);
720 exit(1);
721 }
722 if (!(dseto = SUMA_morphDsetToStd ( dset, M2M,
723 Opt->NearestNode == 1 ? 1:0))) {
724 SUMA_S_Errv("Failed to map %s\n", Opt->in_name);
725 exit(1);
726 }
727 if (0 && !strchr(Opt->out_prefix,'/') ){
728 /* Don't know what made me think that appending the
729 path of the input was a good idea... */
730 s1 = SUMA_append_string(
731 SUMA_FnameGet(Opt->in_name,"pa", SUMAg_CF->cwd),
732 Opt->out_prefix);
733 } else {
734 s1 = SUMA_copy_string(Opt->out_prefix);
735 }
736 s2 = SUMA_RemoveDsetExtension_s(
737 SUMA_FnameGet(Opt->in_name,"l",SUMAg_CF->cwd),
738 SUMA_NO_DSET_FORMAT);
739 uname = SUMA_append_extension(s1,s2);
740 SUMA_free(s1); SUMA_free(s2);
741 oname = SUMA_WriteDset_s (uname, dseto, Opt->oform, 1, 1);
742 if (Opt->NodeDbg>= 0) SUMA_S_Notev("Wrote %s\n", oname);
743 if (oname) SUMA_free(oname); oname=NULL;
744 if (uname) SUMA_free(uname); oname=NULL;
745 if (dseto) SUMA_FreeDset(dseto); dseto = NULL;
746 if (dset) SUMA_FreeDset(dset); dset = NULL;
747 }
748
749 SUMA_LH("Forming the remaining output");
750 outptr = fopen(outname,"w");
751 if (!outptr) {
752 SUMA_SL_Err("Failed to open file for output.\n");
753 exit(1);
754 }
755
756 /* first create the header of the output */
757 SS = SUMA_StringAppend(NULL, NULL);
758 SS = SUMA_StringAppend_va(SS,
759 "#Mapping from nodes on surf 1 (S1) to nodes on surf 2 (S2)\n"
760 "# Surf 1 is labeled %s, idcode:%s\n"
761 "# Surf 2 is labeled %s, idcode:%s\n",
762 SO1->Label, SO1->idcode_str, SO2->Label, SO2->idcode_str);
763 icol = 0;
764 SS = SUMA_StringAppend_va(SS, "#Col. %d:\n"
765 "# S1n (or nj): Index of node on S1\n"
766 , icol);
767 ++icol;
768 if (Opt->NearestNode > 1) {
769 SS = SUMA_StringAppend_va(SS,
770 "#Col. %d..%d:\n"
771 "# S2ne_S1n: Indices of %d nodes on S2 \n"
772 "# that are closest neighbors of nj.\n"
773 "# The first index is that of the node on S2 that is closest \n"
774 "# to nj. If -1 then these values should be ignored because\n"
775 "# in such cases, nj's projection failed.\n"
776 , icol, icol+Opt->NearestNode-1, Opt->NearestNode);
777 icol += Opt->NearestNode;
778 SS = SUMA_StringAppend_va(SS,
779 "#Col. %d..%d:\n"
780 "# S2we_S1n: Weights assigned to nodes on surf 2 (S2) \n"
781 "# that are closest neighbors of nj.\n"
782 , icol, icol+Opt->NearestNode-1, Opt->NearestNode);
783 icol += Opt->NearestNode;
784 } else if (Opt->NearestNode == 1) {
785 SS = SUMA_StringAppend_va(SS,
786 "#Col. %d:\n"
787 "# S2ne_S1n: Index of the node on S2 (label:%s idcode:%s)\n"
788 "# that is the closest neighbor of nj.\n"
789 "# If -1 then this value should be ignored because\n"
790 "# nj's projection failed.\n"
791 , icol, SO2->Label, SO2->idcode_str);
792 ++icol;
793 }
794 if (Opt->NearestTriangle) {
795 SS = SUMA_StringAppend_va(SS,
796 "#Col. %d:\n"
797 "# S2t_S1n: Index of the S2 triangle that hosts node nj on S1.\n"
798 "# In other words, nj's closest projection onto S2 falls on \n"
799 "# triangle S2t_S1n\n"
800 "# If -1 then this value should be ignored because \n"
801 "# nj's projection failed.\n"
802 , icol);
803 ++icol;
804 }
805 if (Opt->ProjectionOnMesh) {
806 SS = SUMA_StringAppend_va(SS,
807 "#Col. %d..%d:\n"
808 "# S2p_S1n: Coordinates of projection of nj onto S2\n"
809 , icol, icol+2);
810 icol += 3;
811 }
812 if (Opt->DistanceToMesh) {
813 SS = SUMA_StringAppend_va(SS,
814 "#Col. %d:\n"
815 "# Closest distance from nj to S2\n"
816 , icol);
817 ++icol;
818 }
819 if (Opt->NearestNodeCoords) {
820 SS = SUMA_StringAppend_va(SS,
821 "#Col. %d .. %d:\n"
822 "# X Y Z coords of nearest node\n"
823 , icol, icol+2);
824 icol += 3;
825 }
826 if (Opt->Data > 0) {
827 if (!Opt->in_name) {
828 SS = SUMA_StringAppend_va(SS,
829 "#Col. %d..%d:\n"
830 "# Interpolation using XYZ coordinates of S2 nodes that neighbor nj\n"
831 "# (same as coordinates of node's projection onto triangle in S2, \n"
832 "# if using barycentric interpolation)\n"
833 , icol, icol+2);
834 icol += 3;
835 } else {
836 SS = SUMA_StringAppend_va(SS,
837 "#Col. %d..%d:\n"
838 "# Interpolation of data at nodes on S2 that neighbor nj\n"
839 "# Data obtained from %s\n"
840 , icol, icol+ncol_data-1, Opt->in_name); icol += ncol_data;
841 }
842 }
843 s = SUMA_HistString("SurfToSurf", argc, argv, NULL);
844 SS = SUMA_StringAppend_va(SS,
845 "#History:\n"
846 "#%s\n", s); SUMA_free(s); s = NULL;
847 SUMA_SS2S(SS,s);
848 fprintf(outptr,"%s\n",s); SUMA_free(s); s = NULL;
849
850 /* put headers atop columns */
851 Nchar = 6; /* if you change this number you'll need to fix formats below */
852 for (i=0; i<icol; ++i) {
853 sprintf(sbuf,"#%s", MV_format_fval2(i, Nchar -1));
854 fprintf(outptr,"%6s ", sbuf);
855 }
856 fprintf(outptr,"\n");
857
858 /* Now put in the values, make sure you parallel columns above! */
859 for (i=0; i<M2M->M1Nn; ++i) {
860 fprintf(outptr,"%6s ", MV_format_fval2(M2M->M1n[i], Nchar));
861 if (Opt->NearestNode > 0) {
862 for (j=0; j<Opt->NearestNode; ++j) {
863 if (j < M2M->M2Nne_M1n[i])
864 fprintf(outptr,"%6s ",
865 MV_format_fval2(M2M->M2ne_M1n[i][j], Nchar));
866 else fprintf(outptr,"%6s ", "-1");
867 } /* Neighboring nodes */
868 }
869 if (Opt->NearestNode > 1) { /* add the weights */
870 for (j=0; j<Opt->NearestNode; ++j) {
871 if (j < M2M->M2Nne_M1n[i])
872 fprintf(outptr,"%6s ",
873 MV_format_fval2(M2M->M2we_M1n[i][j], Nchar));
874 else fprintf(outptr,"%6s ", "0.0");
875 }
876 }
877 if (Opt->NearestTriangle) {
878 fprintf(outptr,"%6s ", MV_format_fval2(M2M->M2t_M1n[i], Nchar));
879 }
880 if (Opt->ProjectionOnMesh) {
881 fprintf(outptr,"%6s ", MV_format_fval2(M2M->M2p_M1n[3*i], Nchar));
882 fprintf(outptr,"%6s ", MV_format_fval2(M2M->M2p_M1n[3*i+1], Nchar));
883 fprintf(outptr,"%6s ", MV_format_fval2(M2M->M2p_M1n[3*i+2], Nchar));
884 }
885 if (Opt->DistanceToMesh) {
886 fprintf(outptr,"%6s ", MV_format_fval2(M2M->PD[i], Nchar));
887 }
888 if (Opt->NearestNodeCoords) {
889 float x=0.0,y=0.0,z=0.0;
890 int n = M2M->M2ne_M1n[i][0];
891 if (n>0) {
892 n = n * SO2->NodeDim;
893 x = SO2->NodeList[n];
894 y = SO2->NodeList[n+1];
895 z = SO2->NodeList[n+2];
896 }
897 fprintf(outptr,"%6s ", MV_format_fval2(x, Nchar));
898 fprintf(outptr,"%6s ", MV_format_fval2(y, Nchar));
899 fprintf(outptr,"%6s ", MV_format_fval2(z, Nchar));
900 }
901 if (dt && Opt->Data > 0) {
902 if (!Opt->in_name) {
903 fprintf(outptr,"%6s ", MV_format_fval2(dt[3*i], Nchar));
904 fprintf(outptr,"%6s ", MV_format_fval2(dt[3*i+1], Nchar));
905 fprintf(outptr,"%6s ", MV_format_fval2(dt[3*i+2], Nchar));
906 } else { /* Column major business */
907 for (j=0; j<ncol_data; ++j) {
908 fprintf(outptr,"%6s ",
909 MV_format_fval2(dt[i+j*M2M->M1Nn], Nchar)); }
910 }
911 }
912 fprintf(outptr,"\n");
913 }
914
915 /* do they want an output surface ? */
916 if (SO_name) {
917 float *tmpfv = NULL;
918 SUMA_LH("Writing surface");
919 tmpfv = SO1->NodeList;
920 SO1->NodeList = dt;
921 if (!SUMA_Save_Surface_Object (SO_name, SO1,
922 ps->o_FT[0], ps->o_FF[0], NULL)) {
923 SUMA_S_Err("Failed to write surface object.\n");
924 exit (1);
925 }
926 SO1->NodeList = tmpfv; tmpfv = NULL;
927 }
928
929 if (N_Spec) {
930 int k=0;
931 for (k=0; k<N_Spec; ++k) {
932 if (!SUMA_FreeSpecFields(&(Spec[k]))) {
933 SUMA_S_Err("Failed to free spec fields");
934 }
935 }
936 SUMA_free(Spec); Spec = NULL; N_Spec = 0;
937 }
938
939 if (projdir) SUMA_free(projdir); projdir = NULL;
940 if (SO_name) SUMA_free(SO_name); SO_name = NULL;
941 if (outptr) fclose(outptr); outptr = NULL;
942 if (dt) SUMA_free(dt); dt = NULL;
943 if (s) SUMA_free(s); s = NULL;
944 if (im_data) mri_free(im_data); im_data = NULL; /* done with the data */
945 if (nodeind) SUMA_free(nodeind); nodeind = NULL;
946 if (M2M) M2M = SUMA_FreeM2M(M2M);
947 if (SO1) SUMA_Free_Surface_Object(SO1); SO1 = NULL;
948 if (SO2) SUMA_Free_Surface_Object(SO2); SO2 = NULL;
949 if (Spec) SUMA_free(Spec); Spec = NULL;
950 if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
951 if (Opt) Opt = SUMA_Free_Generic_Prog_Options_Struct(Opt);
952 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
953 exit(0);
954
955 }
956