1 /* thd_atlas.c */
2 /* functions to access atlas, template, space of datasets */
3
4 /* access current space of dataset */
5
6 #include "mrilib.h"
7 #define SUMA_noFunc
8 #include "suma_objs.h" /* 21 Apr 2020 */
9 /*------------------------------------------------------------*/
10
11 extern int * SUMA_Dijkstra_generic (int N_Node,
12 float *NodeList, int NodeDim, int dist_metric,
13 int *N_Neighbv, int **FirstNeighb, float **FirstNeighbDist,
14 int Nx, int Ny,
15 byte *isNodeInMeshp,
16 int *N_isNodeInMesh, int Method_Number,
17 float *Lfinal, int *N_Path,
18 int verb);
19 static void niml_to_atlas_list(ATLAS_POINT_LIST *atp, char *atlas_file);
20 static void adjust_atlas_point_list(ATLAS_POINT_LIST *atp, char *match_str,
21 float addval);
22 static ATLAS_POINT_LIST *AFNI_atlas_list_to_atlas_point_list(
23 ATLAS_POINT *afni_at_pts, int npts);
24 static int **FirstNeighb=NULL;
25 static float **FirstNeighbDist=NULL;
26 static int *N_Neighb = NULL;
27 static char *jumpspace_name = NULL;
28
THD_get_space(THD_3dim_dataset * dset)29 char * THD_get_space(THD_3dim_dataset *dset)
30 {
31
32 ENTRY("THD_get_space");
33
34 #ifdef DEBUG_SPACES
35 fprintf(stderr,"getting space for dset: %s\n", dset->dblk->diskptr->brick_name);
36 #endif
37
38 if(!ISVALID_DSET(dset)) RETURN(NULL);
39
40 if(dset->atlas_space[0] != '\0'){
41 #ifdef DEBUG_SPACES
42 fprintf(stderr,"Space is %s. Assigned already.\n",dset->atlas_space);
43 #endif
44 RETURN(dset->atlas_space);
45 }
46 switch(dset->view_type){
47 default:
48 case 0:
49 MCW_strncpy(dset->atlas_space, "ORIG", THD_MAX_NAME);
50 break;
51 case 1:
52 MCW_strncpy(dset->atlas_space, "ACPC", THD_MAX_NAME);
53 break;
54 /* view is +tlrc */
55 case 2:
56 /* check environment for default space */
57 /* possible values are TLRC, MNI_ANAT, and MNI */
58 /* but others are possible */
59 MCW_strncpy(dset->atlas_space,
60 TT_whereami_default_spc_name(), THD_MAX_NAME);
61 break;
62 }
63 #ifdef DEBUG_SPACES
64 fprintf(stderr,"Space is %s\n",dset->atlas_space);
65 #endif
66
67 RETURN(dset->atlas_space);
68 }
69
70 /* return the generic space associated with the space of a dataset
71 the principal goal is to give generic TLRC for all flavors of TLRC, TT_N27
72 or generic MNI for all flavors of MNI, MNI_FSL, MNI_ANAT */
THD_get_generic_space(THD_3dim_dataset * dset)73 char * THD_get_generic_space(THD_3dim_dataset *dset)
74 {
75 char *spcstr=NULL, *genspcstr=NULL;
76
77 ENTRY("THD_get_generic_space");
78
79 if(!ISVALID_DSET(dset)) RETURN(NULL);
80 spcstr = THD_get_space(dset); /* space from dataset structure - do not free */
81 if(spcstr!=NULL && *spcstr!='\0')
82 genspcstr = gen_space_str(spcstr); /* space string from space structure - also do not free */
83 if(genspcstr)
84 RETURN(genspcstr);
85 else
86 RETURN(spcstr);
87 }
88
89 /* return the VIEW (ORIG,ACPC or TLRC) for a dataset */
90 /* This is intended as the AFNI output view for a given dataset,
91 i.e.what dataset view should be associated with an AFNI output file.
92 Given an AFNI format input dataset, this will usually be the same
93 as the view of the input (dset->dblk->diskptr->viewcode), but is not
94 guaranteed to be the same. The function solves the problem of
95 NIFTI or other format input and naming the AFNI format output.
96
97 For NIFTI format, the sform or qform codes determine whether the
98 data is in TLRC, MNI or some other aligned space. If it is any aligned
99 space, the dataset will be assigned a TLRC view and a space if one
100 is already defined in an AFNI extension. Otherwise, the TLRC space
101 is assigned to the dataset */
THD_get_view_space(THD_3dim_dataset * dset)102 char * THD_get_view_space(THD_3dim_dataset *dset)
103 {
104 char *spcstr=NULL, *space=NULL;
105
106 ENTRY("THD_get_view_space");
107
108 if(!ISVALID_DSET(dset)) RETURN(NULL);
109 spcstr = dset->dblk->diskptr->viewcode;
110 if(spcstr != NULL)
111 RETURN(spcstr);
112
113 spcstr = THD_get_generic_space(dset); /* space from dataset structure - do not free */
114
115 if (strcmp(spcstr, "ORIG")==0)
116 RETURN("ORIG");
117 if (strcmp(spcstr, "ACPC")==0)
118 RETURN("ACPC");
119 /* all other spaces are assumed to be TLRC */
120 RETURN("TLRC");
121 }
122
123 /* assign space codes used by whereami for specific atlases */
124 int
THD_space_code(char * space)125 THD_space_code(char *space)
126 {
127 ENTRY("THD_space_code");
128
129 if (wami_verb()) {
130 WARNING_message("Better not use codes anymore");
131 }
132 if (strcmp(space, "TLRC")==0)
133 RETURN(AFNI_TLRC_SPC);
134 if (strcmp(space, "MNI_ANAT")==0)
135 RETURN(MNI_ANAT_SPC);
136 if (strcmp(space, "MNI")==0)
137 RETURN(MNI_SPC);
138 RETURN(UNKNOWN_SPC); /* if none of the above, non-standard space */
139 }
140
141
142 /* apply attribute space to dataset structure */
143 void
THD_apply_space_atr(THD_3dim_dataset * dset)144 THD_apply_space_atr(THD_3dim_dataset *dset)
145 {
146 /* allow for non-existent dset structures*/
147 /* allow for environment variable */
148 }
149
150 /* assign space to dataset structure */
151 void
THD_assign_space(THD_3dim_dataset * dset,char * space)152 THD_assign_space(THD_3dim_dataset *dset, char * space)
153 {
154
155 }
156
157 /* make space attribute for dataset */
158 void
THD_make_space_atr(THD_3dim_dataset * dset)159 THD_make_space_atr(THD_3dim_dataset *dset)
160 {
161 /* allow for non-existent dset structures*/
162 /* allow for environment variable */
163
164 }
165
166
167 /* atlas dataset structures */
168 /* atlas datasets will have NIML element attributes that contain
169 segmentation information via index values and segment labels (ROI labels),
170 color scale information that colorizes each segment.
171 This info is along with the space info that will be part of usual datasets. */
172 /* need to also consider the probabilistic atlas type here too */
173 /* read atlas labels */
174
175 /* read atlas color scale */
read_atlas_lut(THD_3dim_dataset * dset)176 ATLAS_LUT *read_atlas_lut(THD_3dim_dataset *dset)
177 {
178 ATLAS_LUT *atlasdset_lut;
179 void *lut_atr;
180
181 ENTRY("read_atlas_lut");
182
183 /* get lookup table from dataset */
184 lut_atr = THD_find_atr(dset->dblk, "ATLAS_LUT");
185 if(lut_atr) {
186 atlasdset_lut = (ATLAS_LUT *)calloc(1,sizeof(ATLAS_LUT));
187 if(atlasdset_lut==NULL) {
188 WARNING_message("Could not allocate memory for Atlas LUT\n");
189 RETURN(NULL);
190 }
191
192 atlasdset_lut->rgblist = calloc(MAXINT, 3);
193 /* memcpy( atlasdset_lut, lut_atr->in , lut_atr->nin ) ;*/
194 RETURN(atlasdset_lut);
195 }
196 RETURN(NULL);
197 }
198
199 /* free list of atlas points */
200 void
free_atlas_point_list(ATLAS_POINT_LIST * apl)201 free_atlas_point_list(ATLAS_POINT_LIST *apl)
202 {
203 if(apl==NULL)
204 return;
205 if(wami_verb() > 1) {
206 INFO_message("Freeing atlas point list with %d points", apl->n_points);
207 print_atlas_point_list(apl);
208 }
209
210 /* if there were any points, free the whole structure too */
211 if(apl->n_points >= 1)
212 free(apl->at_point);
213 free(apl);
214 }
215
216 /* print list of atlas points*/
217 void
print_atlas_point_list(ATLAS_POINT_LIST * apl)218 print_atlas_point_list(ATLAS_POINT_LIST *apl)
219 {
220 int i;
221 ATLAS_POINT *ap;
222 INFO_message("----- Atlas point list: -------");
223
224 if(apl==NULL)
225 return;
226 for(i=0;i<apl->n_points;i++) {
227 ap = apl->at_point+i;
228 INFO_message("%d: \"%s\", \"%s\""
229 " %5.1f %5.1f %5.1f\n",
230 ap->tdval, ap->name, ap->sblabel,
231 ap->xx, ap->yy, ap->zz);
232 }
233 INFO_message("");
234 }
235
236 /* return atlas point from list of atlas points with name of label*/
237 ATLAS_POINT *
atlas_point_named(ATLAS_POINT_LIST * apl,char * name)238 atlas_point_named(ATLAS_POINT_LIST *apl, char *name)
239 {
240 int i;
241 ATLAS_POINT *ap;
242
243 if(apl==NULL)
244 return(NULL);
245 for(i=0;i<apl->n_points;i++) {
246 ap = apl->at_point+i;
247 if(strcmp(ap->name, name)==0)
248 return(ap);
249 }
250 return(NULL);
251 }
252
253 /* return long name for atlas point with name provided */
254 char *
atlas_point_long_name_named(ATLAS_POINT_LIST * apl,char * name)255 atlas_point_long_name_named(ATLAS_POINT_LIST *apl, char *name)
256 {
257 ATLAS_POINT *ap;
258
259 ap = atlas_point_named(apl, name);
260 if(ap) return(ap->longname);
261 return(NULL);
262 }
263
264
265 /* convert a NIML table from a dataset to an atlas list structure */
dset_niml_to_atlas_list(THD_3dim_dataset * dset)266 ATLAS_POINT_LIST * dset_niml_to_atlas_list(THD_3dim_dataset *dset)
267 {
268 ATLAS_POINT_LIST *apl=NULL;
269 int LocalHead = wami_lh();
270 NI_element *nel=NULL;
271 NI_group *ngr=NULL;
272 ATR_string *atr=NULL;
273
274 ENTRY("dset_niml_to_atlas_list");
275
276 if (LocalHead) fprintf(stderr, "assigning NIML attributes to apl.\n");
277
278 atr = THD_find_string_atr( dset->dblk ,
279 "ATLAS_LABEL_TABLE" ) ;
280
281 if (atr) {
282 if (LocalHead) fprintf(stderr, "Label table found in attributes.\n");
283
284 ngr = NI_read_element_fromstring(atr->ch);
285 if ((ngr==NULL) || (ngr->part_num == 0)) {
286 WARNING_message("** WARNING: Poorly formatted ATLAS_LABEL_TABLE\n");
287 if(ngr) NI_free_element(ngr) ;
288 RETURN(NULL);
289 }
290 apl = niml_atlas_label_table_to_atlas_list(ngr);
291 NI_free_element(ngr) ; ngr = NULL;
292 RETURN(apl);
293 }
294 else {
295 if (LocalHead) fprintf(stderr, "Label table NOT found in attributes.\n");
296 RETURN(NULL);
297 }
298
299 RETURN(NULL);
300 }
301
302 /* convert a NIML table to an atlas list structure */
niml_atlas_label_table_to_atlas_list(NI_group * ngr)303 ATLAS_POINT_LIST * niml_atlas_label_table_to_atlas_list(NI_group *ngr)
304 {
305 ATLAS_POINT_LIST *apl;
306 int i, tdlev, tdval;
307 char *temp_str;
308 ATLAS_POINT *at_pt;
309 NI_element *nel=NULL;
310 int LocalHead = wami_lh();
311
312 float cog[3];
313 short okey;
314 ATR_string *atr=NULL;
315
316 ENTRY("niml_atlas_label_table_to_atlas_list");
317
318 if (!ngr) RETURN(NULL);
319
320 /* get each segmented region - the "atlas point" from
321 a NIML formatted string */
322 apl = (ATLAS_POINT_LIST *) calloc(1, sizeof(ATLAS_POINT_LIST));
323 /* assume the number of elements in the group is the number of structures */
324 apl->n_points = ngr->part_num;
325 apl->at_point = (ATLAS_POINT *) calloc(apl->n_points,sizeof(ATLAS_POINT));
326 if(apl->at_point == NULL) {
327 WARNING_message("** WARNING: Poorly formatted ATLAS_POINT_LIST\n");
328 free(apl);
329 RETURN(NULL);
330 }
331
332 for(i=0;i<apl->n_points;i++){
333 /* read the NIML string from the NIML file */
334 nel = (NI_element *)ngr->part[i] ;
335 /* nel = SUMA_FindNgrDataElement(ngr, "ATLAS_POINT", "atlas_point");*/
336 if(!nel) {
337 NI_free_element(ngr) ;
338 RETURN(NULL);
339 }
340 NI_GET_INT(nel, "VAL", tdval);
341 NI_GET_INT(nel, "OKEY", okey);
342 NI_GET_INT(nel, "GYoAR", tdlev);
343 NI_GET_FLOATv(nel, "COG", cog, 3, 0);
344
345 temp_str = NI_get_attribute(nel, "STRUCT");
346 /* update the pointer for the current atlas point */
347 at_pt = &apl->at_point[i];
348
349 /* copy "STRUCT" name - segmentation name */
350 NI_strncpy(at_pt->name,temp_str,ATLAS_CMAX);
351
352 /* mod - drg 09/25/2018 */
353 /* allow for longer name - later add hierarchial name for larger region*/
354 temp_str = NI_get_attribute(nel, "LONGNAME");
355 /* copy "LONGNAME" name - segmentation name */
356 if(temp_str==NULL)
357 NI_strncpy(at_pt->longname,"",ATLAS_CMAX);
358 else
359 NI_strncpy(at_pt->longname,temp_str,ATLAS_CMAX);
360
361 /* sub-brick label for probability maps */
362 temp_str = NI_get_attribute(nel, "SB_LABEL");
363 if(temp_str==NULL)
364 NI_strncpy(at_pt->sblabel,"",ATLAS_CMAX);
365 else
366 NI_strncpy(at_pt->sblabel,temp_str,ATLAS_CMAX);
367
368 at_pt->tdval = tdval;
369 at_pt->okey = okey;
370 at_pt->tdlev = tdlev;
371 at_pt->xx = cog[0];
372 at_pt->yy = cog[1];
373 at_pt->zz = cog[2];
374 }
375
376 RETURN(apl);
377 }
378
379
380 /* search niml atlas label attribute for label
381 * if label matches regular or longname, return index of structure
382 * atlas does not need to be in atlas list */
dset_atlas_label_index(THD_3dim_dataset * dset,char * label)383 int dset_atlas_label_index(THD_3dim_dataset *dset, char *label)
384 {
385 ATLAS_POINT_LIST *apl=NULL;
386 int LocalHead = wami_lh();
387 NI_element *nel=NULL;
388 NI_group *ngr=NULL;
389 ATR_string *atr=NULL;
390
391 int i, tdlev, tdval, nr;
392 char *temp_str1, *temp_str2;
393 ATLAS_POINT *at_pt;
394
395 ENTRY("dset_atlas_label_index");
396
397 atr = THD_find_string_atr( dset->dblk ,
398 "ATLAS_LABEL_TABLE" ) ;
399
400 if (atr) {
401 if (LocalHead) fprintf(stderr, "Label table found in attributes.\n");
402
403 ngr = NI_read_element_fromstring(atr->ch);
404 if ((ngr==NULL) || (ngr->part_num == 0)) {
405 WARNING_message("** WARNING: Poorly formatted ATLAS_LABEL_TABLE\n");
406 if(ngr) NI_free_element(ngr) ;
407 RETURN(0);
408 }
409 }
410 else {
411 RETURN(0);
412 }
413
414 if (!ngr) RETURN(0);
415
416 if(LocalHead)
417 printf("Could read atlas labels in attributes\n");
418
419 /* get each segmented region - the "atlas point" from
420 a NIML formatted string */
421 apl = (ATLAS_POINT_LIST *) calloc(1, sizeof(ATLAS_POINT_LIST));
422 /* assume the number of elements in the group is the number of structures */
423 apl->n_points = ngr->part_num;
424 apl->at_point = (ATLAS_POINT *) calloc(apl->n_points,sizeof(ATLAS_POINT));
425 if(apl->at_point == NULL) {
426 free(apl);
427 RETURN(0);
428 }
429 if(LocalHead)
430 printf("Number of points in atlas is %d\n", apl->n_points);
431
432 /* replace spaces in label with underscores */
433 nr = str_replace_char(label, ' ', '_');
434
435 for(i=0;i<apl->n_points;i++){
436 /* read the NIML string from the NIML file */
437 nel = (NI_element *)ngr->part[i] ;
438 if(!nel) {
439 free(apl);
440 NI_free_element(ngr) ;
441 RETURN(0);
442 }
443 NI_GET_INT(nel, "VAL", tdval);
444
445 temp_str1 = NI_get_attribute(nel, "STRUCT");
446 temp_str2 = NI_get_attribute(nel, "LONGNAME");
447 if(LocalHead)
448 printf("index %d STRUCT: %s LONGNAME %s\n", tdval, temp_str1, temp_str2);
449 // (if(temp_str2)(strcmp(label, temp_str2)==0))
450 /* match either short or longname */
451 /* replace spaces in label with underscores-easier to pass in scripts */
452 if(temp_str1) nr = str_replace_char(temp_str1, ' ', '_');
453 if(temp_str1 && strcmp(label, temp_str1)==0) {
454 if(LocalHead)
455 printf("Found label %s at index %d\n", label, tdval);
456
457 free(apl);
458 NI_free_element(ngr) ;
459 RETURN(tdval);
460 }
461 else {
462 if(temp_str2) nr = str_replace_char(temp_str1, ' ', '_');
463 if(temp_str2 && strcmp(label, temp_str2)==0) {
464 if(LocalHead)
465 printf("Found label %s at index %d.\n"
466 "Replaced %d spaces with underscores\n", label, tdval,nr);
467
468 free(apl);
469 NI_free_element(ngr) ;
470 RETURN(tdval);
471 }
472 }
473 }
474 free(apl);
475 NI_free_element(ngr);
476
477 RETURN(0);
478 }
479
identity_xform_chain(char * space_name)480 ATLAS_XFORM_LIST *identity_xform_chain(char *space_name)
481 {
482 ATLAS_XFORM *xf;
483 ATLAS_XFORM_LIST *xfl=NULL;
484
485 xf = identity_xform(); /* assign identity matrix */
486 free(xf->source); free(xf->dest);
487 xf->source = nifti_strdup(space_name);
488 xf->dest = nifti_strdup(space_name);
489 xfl = (ATLAS_XFORM_LIST *) calloc(1, sizeof(ATLAS_XFORM));
490 xfl->xform = xf;
491 xfl->nxforms = 1;
492 return(xfl);
493 }
494
495
496 /* write atlas labels */
497
498 /* write atlas color scale */
499
500 /* read transformation for space 1 to space 2 */
501
502 /* write transformation for space 1 to space 2 */
503
504
505 /* report xform chain to go between src and dest spaces */
report_xform_chain(char * src,char * dest,int report)506 ATLAS_XFORM_LIST *report_xform_chain(char *src, char *dest, int report)
507 {
508 ATLAS_XFORM_LIST *xfl;
509 int srci, desti;
510 ATLAS_SPACE_LIST *asl=get_G_space_list();
511
512 if (!src || !dest) return(NULL);
513
514 /* atlas testing */
515 if (!strcmp(src, dest)) {
516 if(wami_verb() > 1)
517 INFO_message("Chain is from and to same space %s", src);
518 xfl = identity_xform_chain(src);
519 } else {
520 srci = find_atlas_space_index(src);
521 if(srci<0 && wami_verb()){
522 INFO_message("Could not find source space %s in database", src);
523 }
524 desti = find_atlas_space_index(dest);
525 if(desti<0 && wami_verb()){
526 INFO_message("Could not find destination space %s in database", dest);
527 print_space_list(NULL);
528 }
529 /* check if we're going nowhere */
530 if(srci==desti) {
531 if(wami_verb() > 1)
532 INFO_message("Chain is from and to same space");
533 xfl = identity_xform_chain(src);
534 } else {
535 if (srci >= 0 && desti >= 0) {
536 xfl = get_xform_chain(asl->space+srci,
537 asl->space+desti);
538 } else {
539 xfl = NULL;
540 }
541 }
542 }
543 if (report) print_xform_list(xfl);
544 return(xfl);
545 }
546
547 /* report which spaces are available to go between src and all other spaces */
report_available_spaces(char * src)548 void report_available_spaces(char *src)
549 {
550 ATLAS_SPACE_LIST *spl;
551
552 spl = find_available_spaces(src, NULL);
553 print_space_list(spl);
554 free_space_list(spl);
555 }
556
557
558 /* initialize space structures for dealing with templates and atlases */
init_space_structs(ATLAS_XFORM_LIST ** atlas_xfl,ATLAS_LIST ** atlas_alist,ATLAS_SPACE_LIST ** atlas_spaces,ATLAS_TEMPLATE_LIST ** atlas_templates)559 int init_space_structs(ATLAS_XFORM_LIST **atlas_xfl,
560 ATLAS_LIST **atlas_alist,
561 ATLAS_SPACE_LIST **atlas_spaces,
562 ATLAS_TEMPLATE_LIST **atlas_templates)
563 {
564 *atlas_alist = (ATLAS_LIST *) calloc(1,sizeof(ATLAS_LIST));
565 *atlas_spaces = (ATLAS_SPACE_LIST *) calloc(1,sizeof(ATLAS_SPACE_LIST));
566 *atlas_templates = (ATLAS_TEMPLATE_LIST *)
567 calloc(1, sizeof(ATLAS_TEMPLATE_LIST));
568 *atlas_xfl = (ATLAS_XFORM_LIST *) calloc(1, sizeof(ATLAS_XFORM_LIST));
569 (*atlas_xfl)->nxforms = 0;
570 (*atlas_xfl)->xform = NULL;
571 (*atlas_alist)->natlases = 0;
572 (*atlas_spaces)->nspaces = 0;
573 (*atlas_templates)->ntemplates = 0;
574
575 if( !(*atlas_alist) || !(*atlas_spaces) ||
576 !(*atlas_templates) || !(*atlas_xfl)) {
577 return(0);
578 }
579 return(1);
580 }
581
582
583 /* read niml file elements into C structures */
read_space_niml_file(char * fname,ATLAS_XFORM_LIST * atlas_xfl,ATLAS_LIST * atlas_alist,ATLAS_SPACE_LIST * atlas_spaces,ATLAS_TEMPLATE_LIST * atlas_templates,THD_string_array * sar)584 int read_space_niml_file( char *fname, ATLAS_XFORM_LIST *atlas_xfl,
585 ATLAS_LIST *atlas_alist,
586 ATLAS_SPACE_LIST *atlas_spaces,
587 ATLAS_TEMPLATE_LIST *atlas_templates,
588 THD_string_array *sar)
589 {
590 NI_element *nel;
591 NI_stream space_niml;
592 char *fnamet=NULL, *fnameclean=NULL;
593 int found = 0;
594
595 if (!fname) return(found);
596
597 if (!(fnameclean = af_strnstr(fname, "file:", 5))) {
598 fnamet = (char *)calloc(6+strlen(fname), sizeof(char));
599 sprintf(fnamet,"file:%s", fname);
600 space_niml = NI_stream_open(fnamet,"r");
601 free(fnamet);
602 fnameclean = fname;
603 } else {
604 space_niml = NI_stream_open(fname,"r");
605 }
606
607 if (!space_niml) {
608 ERROR_message("Failed to NI_stream open %s\n", fnameclean);
609 return(found);
610 }
611
612 nel = (NI_element *) 1;
613 while(nel) {
614 if(wami_verb() > 2)
615 INFO_message("reading elements\n");
616 if ((nel = NI_read_element(space_niml, 100))) {
617 found += add_atlas_nel(nel, atlas_xfl,
618 atlas_alist,atlas_spaces,atlas_templates, sar,
619 THD_filepath(fnameclean));
620 NI_free_element(nel); /* don't need the NIML element anymore */
621 }
622 }
623
624 NI_stream_close(space_niml);
625
626 return(found);
627 }
628
add_atlas_nel(NI_element * nel,ATLAS_XFORM_LIST * atlas_xfl,ATLAS_LIST * atlas_alist,ATLAS_SPACE_LIST * atlas_spaces,ATLAS_TEMPLATE_LIST * atlas_templates,THD_string_array * sar,char * parentdir)629 int add_atlas_nel(NI_element *nel, ATLAS_XFORM_LIST *atlas_xfl,
630 ATLAS_LIST *atlas_alist,
631 ATLAS_SPACE_LIST *atlas_spaces,
632 ATLAS_TEMPLATE_LIST *atlas_templates,
633 THD_string_array *sar,
634 char *parentdir) {
635
636 int found = 0;
637
638 if (!nel) return(found);
639
640 if(wami_verb() > 2)
641 INFO_message("nel name %s\n", nel->name);
642 if (nel->type == NI_ELEMENT_TYPE) {
643 if(strcmp(nel->name, "TEMPLATE_SPACE") == 0) {
644 atlas_spaces->nspaces++;
645 if(wami_verb() > 1){
646 INFO_message("Template space\n"
647 "number of spaces now %d\n",
648 atlas_spaces->nspaces);
649 }
650 if(atlas_spaces->nspaces==1) {
651 if(wami_verb() > 2)
652 INFO_message("initial memory allocation for spaces");
653 atlas_spaces->space =
654 (ATLAS_SPACE *) calloc(1, sizeof(ATLAS_SPACE));
655 } else {
656 atlas_spaces->space =
657 (ATLAS_SPACE *) realloc( atlas_spaces->space,
658 atlas_spaces->nspaces * sizeof(ATLAS_SPACE));
659 }
660 atlas_read_atlas_space(
661 nel, &atlas_spaces->space[atlas_spaces->nspaces-1]);
662
663 found = 1;
664 }
665 if(strcmp(nel->name, "XFORM") == 0) {
666 atlas_xfl->nxforms++;
667 if(wami_verb() > 2){
668 INFO_message("space XFORM\n");
669 INFO_message("number of xforms now %d\n", atlas_xfl->nxforms);
670 }
671 if(atlas_xfl->nxforms==1){
672 if(wami_verb() > 2)
673 INFO_message("initial memory allocation for xforms");
674 atlas_xfl->xform =
675 (ATLAS_XFORM *) calloc(1, sizeof(ATLAS_XFORM));
676 }
677 else
678 atlas_xfl->xform = (ATLAS_XFORM *) realloc(
679 atlas_xfl->xform,
680 atlas_xfl->nxforms * sizeof(ATLAS_XFORM));
681 atlas_read_xform(nel, &atlas_xfl->xform[atlas_xfl->nxforms-1]);
682 found = 1;
683 }
684 if(strcmp(nel->name, "ATLAS") == 0) {
685 atlas_alist->natlases++;
686 if(wami_verb() > 2){
687 INFO_message("Number of atlases now %d\n",
688 atlas_alist->natlases);
689 }
690 if(atlas_alist->natlases==1){
691 if(wami_verb() > 2)
692 INFO_message("initial memory allocation for atlases");
693 atlas_alist->atlas = (ATLAS *) calloc(1, sizeof(ATLAS));
694 } else {
695 atlas_alist->atlas = (ATLAS *) realloc(
696 atlas_alist->atlas, atlas_alist->natlases * sizeof(ATLAS));
697 memset(&(atlas_alist->atlas[atlas_alist->natlases-1]), 0,
698 sizeof(ATLAS));
699 }
700 atlas_read_atlas(nel,
701 &atlas_alist->atlas[atlas_alist->natlases-1],
702 parentdir);
703 if (sar &&
704 atlas_alist->atlas[atlas_alist->natlases-1].name){
705 ADDUTO_SARR(sar,
706 atlas_alist->atlas[atlas_alist->natlases-1].name);
707 } /* add the name of that atlas to the string array */
708 found = 1;
709 }
710
711 if(strcmp(nel->name, "TEMPLATE") == 0) {
712 atlas_templates->ntemplates++;
713 if(wami_verb() > 2){
714 INFO_message("Atlas template\n");
715 INFO_message("number of templates now %d\n",
716 atlas_templates->ntemplates);
717 }
718 if(atlas_templates->ntemplates==1){
719 if(wami_verb() > 2)
720 INFO_message("initial memory allocation for templates");
721 atlas_templates->atlas_template =
722 (ATLAS_TEMPLATE *) calloc(1,sizeof(ATLAS_TEMPLATE));
723 }
724 else
725 atlas_templates->atlas_template = (ATLAS_TEMPLATE *) realloc(
726 atlas_templates->atlas_template,
727 atlas_templates->ntemplates * sizeof(ATLAS_TEMPLATE));
728 atlas_read_template(nel,
729 &atlas_templates->atlas_template[
730 atlas_templates->ntemplates-1]);
731 found = 1;
732 }
733 }
734 return(found);
735 }
736
737 /* compute what spaces are neighbors of the others */
738 /* in preparation for Dijkstra search */
make_space_neighborhood(ATLAS_SPACE_LIST * at_spl,ATLAS_XFORM_LIST * atlas_xfl)739 int make_space_neighborhood( ATLAS_SPACE_LIST *at_spl,
740 ATLAS_XFORM_LIST *atlas_xfl)
741 {
742 /* See Ziad's NIH-5 labbook pp 46 for graph */
743
744
745 int i, j, nspaces, inv_xf, neighbor_i;
746 ATLAS_SPACE *atlas_space, *dest_space;
747 ATLAS_XFORM *xform;
748
749 nspaces = at_spl->nspaces;
750 if(nspaces == 0) {
751 if(wami_verb() > 1)
752 INFO_message("no spaces to compute paths among");
753 FirstNeighb = NULL; FirstNeighbDist = NULL; N_Neighb = NULL;
754 return(-1);
755 }
756 FirstNeighb = (int **) calloc(nspaces , sizeof(int *));
757 FirstNeighbDist = (float **) calloc(nspaces , sizeof(float *));
758 N_Neighb = (int *) calloc(nspaces , sizeof(int));
759
760 if(wami_verb() > 2)
761 INFO_message("initial memory allocation for neighbors: nspaces = %d",
762 nspaces);
763 if((FirstNeighb==NULL) || (FirstNeighbDist==NULL) ||
764 (N_Neighb==NULL)) {
765 WARNING_message("Could not allocate space for atlas neighborhood.");
766 return(-1);
767 }
768 /* loop through all the spaces and see if they can be directly transformed
769 (or inversely transformed) to any other space in the list of spaces */
770 for(i=0;i<nspaces;i++){
771 neighbor_i = 0;
772 atlas_space = at_spl->space+i;
773 for(j=0;j<nspaces;j++) {
774 dest_space = at_spl->space+j;
775 if(wami_verb() > 1)
776 INFO_message("Computing path from %s(%d) to %s(%d)",
777 atlas_space->atlas_space, i,
778 dest_space->atlas_space, j);
779
780 if(i==j) continue; /* don't need to match the same space */
781 xform = get_xform_neighbor(
782 atlas_xfl, atlas_space, dest_space, &inv_xf);
783 if(xform!=NULL){
784 if(neighbor_i==0){
785 FirstNeighb[i] = (int *) calloc(1, sizeof(int));
786 FirstNeighbDist[i] = (float *) calloc(1, sizeof(float));
787 }
788 else {
789 FirstNeighb[i] = (int *) realloc(FirstNeighb[i],
790 (neighbor_i+1)*sizeof(int));
791 FirstNeighbDist[i] = (float *) realloc(FirstNeighbDist[i],
792 (neighbor_i+1)*sizeof(float));
793 }
794 if((FirstNeighb[i]==NULL) || (FirstNeighbDist[i]==NULL)) {
795 WARNING_message("Could not allocate space for "
796 "atlas neighborhood");
797 return(-1);
798 }
799
800 FirstNeighb[i][neighbor_i] = j;
801 FirstNeighbDist[i][neighbor_i++] = xform->dist;
802 if(wami_verb() > 1){
803 INFO_message("neighbor found for space %d with space %d",i,j);
804 INFO_message("xform %s with dist %f",
805 xform->xform_name, xform->dist);
806 }
807 }
808 }
809 N_Neighb[i] = neighbor_i;
810 }
811
812 return(0);
813 }
814
815 /* search through list of transformations for direct or inverse transformation
816 of source template space to destination template space */
817 ATLAS_XFORM *
get_xform_neighbor(ATLAS_XFORM_LIST * atlas_xfl,ATLAS_SPACE * at_space,ATLAS_SPACE * dest_space,int * inv_xf)818 get_xform_neighbor(ATLAS_XFORM_LIST *atlas_xfl, ATLAS_SPACE *at_space,
819 ATLAS_SPACE *dest_space, int *inv_xf)
820 {
821 int i, cc;
822 char *srcstr, *deststr, *xfsrc, *xfdest;
823 ATLAS_XFORM *xf, *xf2 = NULL;
824
825 srcstr = at_space->atlas_space;
826 deststr = dest_space->atlas_space;
827
828 *inv_xf = 0;
829
830 /* check if xform from src to dest space is in list of xforms */
831 for(i=0;i<atlas_xfl->nxforms;i++) {
832 xf = atlas_xfl->xform+i;
833 xfsrc = xf->source;
834 xfdest = xf->dest;
835 if((strcmp(srcstr, xfsrc)==0) && (strcmp(deststr,xfdest)==0)) {
836 return(xf);
837 }
838 }
839
840 /* if we made it here, then check for dest to src, and mark to invert */
841 for(i=0;i<atlas_xfl->nxforms;i++) {
842 xf = atlas_xfl->xform+i;
843 xfsrc = xf->source;
844 xfdest = xf->dest;
845 /* check if inverse direction is available */
846 if((strcmp(deststr, xfsrc)==0) && (strcmp(srcstr,xfdest)==0)) {
847 /* is the original xform invertible */
848 xf2 = (ATLAS_XFORM *) calloc(1, sizeof(ATLAS_XFORM));
849 if(copy_xform(xf, xf2)!=0){
850 WARNING_message("Could not create copy of xform for path");
851 return(NULL);
852 }
853 xf2->inverse = 1;
854 cc = invert_xform(xf2);
855 free_xform(xf2);
856 free(xf2);
857 if(cc) {
858 if(wami_verb() > 1){
859 INFO_message("Can not invert transform in path from %s to %s",
860 xfsrc, xfdest);
861 }
862 }
863 else {
864 if(wami_verb() > 1){
865 INFO_message("Using invertible transform in path from %s to %s",
866 xfsrc, xfdest);
867 }
868 *inv_xf = 1;
869 return(xf);
870 }
871 }
872 }
873
874 return(NULL);
875 }
876
877 /* find shortest path to go from one space to a destination space */
878 /* return the list of transformations needed to accomplish this */
get_xform_chain(ATLAS_SPACE * at_space,ATLAS_SPACE * dest_space)879 ATLAS_XFORM_LIST * get_xform_chain( ATLAS_SPACE *at_space,
880 ATLAS_SPACE *dest_space)
881 {
882 int srci, desti;
883 int N_n, kk, *nPath;
884 float nDistance;
885 ATLAS_XFORM_LIST *xfl=NULL;
886 ATLAS_SPACE_LIST *asl=get_G_space_list();
887 ATLAS_XFORM_LIST *axl=get_G_xform_list();
888 /* you might want to return identity right here if
889 at_space->atlas_space == dest_space->atlas_space,
890 even if find_atlas_space can't find them ... */
891
892 /* find index for input spaces */
893 if ((srci = find_atlas_space(asl, at_space))<0) {
894 ERROR_message("input space %s/%s not in atlas_spaces",
895 at_space->atlas_space, at_space->generic_space);
896 print_space_list(asl);
897 return(NULL);
898 }
899 if ((desti = find_atlas_space(asl, dest_space))<0) {
900 ERROR_message("destination space %s/%s not in atlas_spaces",
901 dest_space->atlas_space, dest_space->generic_space);
902 return(NULL);
903 };
904 /* if src and dest are the same, should return identity right away */
905 /* check if neighborhood is defined */
906 if((N_Neighb==NULL) || (FirstNeighbDist==NULL)) return (NULL);
907 /* search neighborhood for shortest path between indices */
908 if ( !(nPath = SUMA_Dijkstra_generic (
909 asl->nspaces,
910 NULL, -1, 0,
911 N_Neighb, FirstNeighb, FirstNeighbDist,
912 srci, desti,
913 NULL, NULL,
914 1,
915 &nDistance, &N_n, 0)) ) {
916 if(wami_verb()>1) fprintf(stderr,
917 "No path found in Dijkstra from %d to %d space", srci, desti);
918 return(NULL); /* no path found */
919 } else {
920 if(wami_verb() > 1){
921 INFO_message("Number of spaces to traverse %d with distance %.2f ",
922 N_n, nDistance);
923 fprintf(stderr, "spaces in chain by index: ");
924 for(kk=0; kk<N_n; ++kk)
925 fprintf(stderr, "%d ", nPath[kk]);
926 fprintf(stderr, "\n");
927 }
928
929 xfl = pathlist_to_xform_list(nPath, N_n,
930 axl, asl);
931 free(nPath); nPath = NULL;
932 }
933
934 return(xfl);
935 }
936
937 /* find space index that matches space name */
find_atlas_space(ATLAS_SPACE_LIST * at_spl,ATLAS_SPACE * at_space)938 int find_atlas_space(ATLAS_SPACE_LIST *at_spl, ATLAS_SPACE *at_space)
939 {
940 int i;
941 ATLAS_SPACE *sp1;
942
943 for(i=0;i<at_spl->nspaces;i++) {
944 sp1 = at_spl->space+i;
945 if(strcmp(sp1->atlas_space, at_space->atlas_space)==0)
946 return(i);
947 }
948 return(-1); /* not found */
949 }
950
951 /* find space index in space list that matches string */
find_atlas_space_index(char * spacename)952 int find_atlas_space_index(char *spacename)
953 {
954 int i;
955 ATLAS_SPACE *sp1;
956 ATLAS_SPACE_LIST *asl=get_G_space_list();
957
958 if((spacename==NULL) || (strcmp(spacename, "")==0)
959 || !asl) {
960 if (wami_verb()) {
961 ERROR_message("Null input: spacename = %s, asl = %p",
962 STR_PRINT(spacename), asl);
963 }
964 return(-1);
965 }
966 for(i=0;i<asl->nspaces;i++) {
967 sp1 = asl->space+i;
968 if(strcmp(sp1->atlas_space, spacename)==0)
969 return(i);
970 }
971 return(-1); /* not found */
972 }
973
974 /* return the atlas_space structure corresponding to the atlas_space name
975 in the dataset structure */
dset_space(THD_3dim_dataset * dset)976 ATLAS_SPACE * dset_space(THD_3dim_dataset *dset)
977 {
978 int spacei;
979 ATLAS_SPACE_LIST *asl=get_G_space_list();
980
981 spacei = find_atlas_space_index(dset->atlas_space);
982 if(spacei==-1)
983 return(NULL);
984
985 return(asl->space+spacei);
986 }
987
988 /* find all available spaces that can be transformed to from some original space */
find_available_spaces(char * src_space_name,ATLAS_SPACE_LIST * this_spl)989 ATLAS_SPACE_LIST *find_available_spaces( char *src_space_name,
990 ATLAS_SPACE_LIST *this_spl)
991 {
992 int i, curr_i, nspaces=0;
993 ATLAS_SPACE_LIST *spl;
994 ATLAS_SPACE_LIST *search_list=NULL;
995 ATLAS_XFORM_LIST *xfl;
996 ATLAS_SPACE *xs, *spl_space, *src_space;
997
998 if (!this_spl) search_list = get_G_space_list();
999 else search_list = this_spl;
1000
1001 spl = NULL;
1002
1003 /* find index of current space string in static space list */
1004 curr_i = find_atlas_space_index(src_space_name);
1005 src_space = search_list->space+curr_i;
1006
1007 /* search through all spaces */
1008 for(i=0;i<search_list->nspaces;i++){
1009 if(i==curr_i) continue; /* don't count the current space */
1010 xs = search_list->space+i; /* pointer to indexed space */
1011 xfl = get_xform_chain(src_space, xs);
1012 if(xfl){ /* found a transformation */
1013 if(wami_verb() > 1)
1014 INFO_message("Found an available space: %s", xs->atlas_space);
1015 free_xform_list(xfl); /* don't actually need the xform list */
1016
1017 if(spl==NULL) {
1018 spl = (ATLAS_SPACE_LIST *) calloc(1, sizeof(ATLAS_SPACE_LIST));
1019 spl->space = (ATLAS_SPACE *) calloc(1, sizeof(ATLAS_SPACE));
1020 nspaces = 1;
1021 }
1022 else {
1023 nspaces++;
1024 spl->space = (ATLAS_SPACE *) realloc(
1025 spl->space, nspaces * sizeof(ATLAS_SPACE));
1026 }
1027
1028 if((spl==NULL)||(spl->space==NULL)) {
1029 WARNING_message("Could not allocate available space transformation");
1030 return(NULL);
1031 }
1032 spl_space = spl->space+nspaces-1;
1033 spl_space->atlas_space = nifti_strdup(xs->atlas_space);
1034 spl_space->generic_space = nifti_strdup(xs->generic_space);
1035
1036 if( (spl_space->atlas_space == NULL) ||
1037 (spl_space->generic_space == NULL)) {
1038 WARNING_message("Could not allocate template space strings");
1039 return(NULL);
1040 }
1041 spl->nspaces = nspaces;
1042
1043 }
1044 }
1045 if(spl) {
1046 spl->nspaces = nspaces;
1047 if(wami_verb() > 1)
1048 INFO_message("There are %d spaces available", spl->nspaces);
1049 } else {
1050 if(wami_verb() > 1) {
1051 print_space_list(search_list);
1052 INFO_message("no spaces reachable from source space: %s",
1053 src_space_name);
1054 }
1055 }
1056
1057 return(spl);
1058 }
1059
pathlist_to_xform_list(int * nPath,int N_n,ATLAS_XFORM_LIST * atlas_xfl,ATLAS_SPACE_LIST * at_spl)1060 ATLAS_XFORM_LIST * pathlist_to_xform_list(int *nPath, int N_n,
1061 ATLAS_XFORM_LIST *atlas_xfl,
1062 ATLAS_SPACE_LIST *at_spl)
1063 {
1064 int kk, inv_xf;
1065 ATLAS_XFORM_LIST *xflc = NULL;
1066 ATLAS_XFORM *a_xform, *xxform;
1067 ATLAS_SPACE *sp1, *sp2;
1068
1069 xflc = (ATLAS_XFORM_LIST *) calloc(1, sizeof(ATLAS_XFORM_LIST));
1070 xflc->xform = (ATLAS_XFORM *) calloc((N_n-1), sizeof(ATLAS_XFORM));
1071 xflc->nxforms = N_n-1;
1072
1073 for(kk=0;kk<N_n-1;++kk) {
1074 sp1 = at_spl->space+nPath[kk]; /* starting space */
1075 sp2 = at_spl->space+nPath[kk+1]; /* next space */
1076 a_xform = get_xform_neighbor(atlas_xfl,sp1, sp2, &inv_xf);
1077 if(wami_verb() > 1){
1078 INFO_message("space%d %s to space%d %s using xform %s",
1079 kk, sp1->atlas_space, kk+1, sp2->atlas_space, a_xform->xform_name);
1080 }
1081
1082 xxform = xflc->xform+kk;
1083 if(copy_xform(a_xform, xxform)!=0){
1084 WARNING_message("Could not create copy of xform for path");
1085 return(NULL);
1086 }
1087
1088
1089 if(inv_xf) /* if inverting xform, take inverse of inverse field */
1090 xxform->inverse = !(a_xform->inverse);
1091
1092 if(wami_verb() > 1){
1093 print_xform(xxform);
1094 }
1095
1096 }
1097
1098 if(wami_verb() > 1){
1099 INFO_message("Number of xforms in chain is %d", xflc->nxforms);
1100 }
1101
1102 return(xflc);
1103 }
1104
1105
1106 /* copy elements of xform structure, allocating space for each element */
1107 /* dest_xform already exists as a structure, it just needs to be populated */
1108 int
copy_xform(ATLAS_XFORM * src_xform,ATLAS_XFORM * dest_xform)1109 copy_xform(ATLAS_XFORM *src_xform, ATLAS_XFORM *dest_xform)
1110 {
1111 memset(dest_xform, 0, sizeof(ATLAS_XFORM));
1112 dest_xform->xform_type = nifti_strdup(src_xform->xform_type);
1113 dest_xform->xform_name = nifti_strdup(src_xform->xform_name);
1114 dest_xform->source = nifti_strdup(src_xform->source);
1115 dest_xform->dest = nifti_strdup(src_xform->dest);
1116 dest_xform->coord_order = nifti_strdup(src_xform->coord_order);
1117
1118 if((dest_xform->xform_type==NULL) ||(dest_xform->xform_name==NULL) ||
1119 (dest_xform->source==NULL) || (dest_xform->dest==NULL) ||
1120 (dest_xform->coord_order==NULL))
1121 return(1);
1122 dest_xform->dist = src_xform->dist;
1123 dest_xform->inverse = src_xform->inverse;
1124 dest_xform->post = src_xform->post;
1125 dest_xform->nelts = src_xform->nelts;
1126 if(dest_xform->nelts==0) return(0);
1127 dest_xform->xform = calloc(dest_xform->nelts, sizeof(float));
1128 if(dest_xform->xform==NULL) return(1);
1129 memcpy(dest_xform->xform, src_xform->xform, dest_xform->nelts*sizeof(float));
1130 return(0);
1131 }
1132
1133 /* copy elements of xform structure, allocating space for each element */
1134 /* dest_xform already exists as a structure, it just needs to be populated */
identity_xform()1135 ATLAS_XFORM *identity_xform()
1136 {
1137 ATLAS_XFORM *dest_xform;
1138 float *fptr;
1139
1140 dest_xform = (ATLAS_XFORM *) calloc(1,sizeof(ATLAS_XFORM));
1141 dest_xform->xform_type = nifti_strdup("Identity");
1142 dest_xform->xform_name = nifti_strdup("Identity");
1143 dest_xform->source = nifti_strdup("");
1144 dest_xform->dest = nifti_strdup("");
1145 dest_xform->coord_order = nifti_strdup("rai");
1146
1147 if((dest_xform->xform_type==NULL) ||(dest_xform->xform_name==NULL) ||
1148 (dest_xform->source==NULL) || (dest_xform->dest==NULL) ||
1149 (dest_xform->coord_order==NULL))
1150 return(NULL);
1151 dest_xform->dist = 0.01;
1152 dest_xform->inverse = 0;
1153 dest_xform->post = 1;
1154
1155 dest_xform->nelts = 1;
1156 if(dest_xform->nelts==0) return(dest_xform);
1157 dest_xform->xform = calloc(dest_xform->nelts, sizeof(float));
1158 if(dest_xform->xform==NULL) return(NULL);
1159 fptr = (float *) dest_xform->xform;
1160 *fptr = 1.0;
1161 return(dest_xform);
1162 }
1163
1164
1165 /* free list of xforms */
free_xform_list(ATLAS_XFORM_LIST * xfl)1166 void free_xform_list(ATLAS_XFORM_LIST *xfl)
1167 {
1168 int i;
1169
1170 if(xfl==NULL)
1171 return;
1172 for(i=(xfl->nxforms)-1;i>=0;i--){
1173 free_xform(xfl->xform+i);
1174 }
1175 free(xfl->xform);
1176 free(xfl);
1177 }
1178
1179 /* free an atlas_xform structure */
free_xform(ATLAS_XFORM * xf)1180 void free_xform(ATLAS_XFORM *xf)
1181 {
1182 if(xf==NULL)
1183 return;
1184 free(xf->xform);
1185 free(xf->xform_type); free(xf->xform_name); free(xf->source); free(xf->dest);
1186 free(xf->coord_order);
1187 }
1188
1189 /* free list of spaces */
free_space_list(ATLAS_SPACE_LIST * xsl)1190 void free_space_list(ATLAS_SPACE_LIST *xsl)
1191 {
1192 int i;
1193
1194 if(xsl==NULL)
1195 return;
1196 for(i=0;i<xsl->nspaces;i++)
1197 free_space(xsl->space+i);
1198 free(xsl->space);
1199 free(xsl);
1200 }
1201
1202 /* free an atlas_space structure */
free_space(ATLAS_SPACE * xs)1203 void free_space(ATLAS_SPACE *xs)
1204 {
1205 if(xs==NULL)
1206 return;
1207 free(xs->atlas_space); free(xs->generic_space);
1208 }
1209
1210 /* free list of spaces */
free_template_list(ATLAS_TEMPLATE_LIST * xtl)1211 void free_template_list(ATLAS_TEMPLATE_LIST *xtl)
1212 {
1213 int i;
1214
1215 if(xtl==NULL)
1216 return;
1217 for(i=0;i<xtl->ntemplates;i++)
1218 free_template(xtl->atlas_template+i);
1219
1220 if(xtl->ntemplates >= 1)
1221 free(xtl->atlas_template);
1222 free(xtl);
1223 }
1224
1225 /* free an atlas_template structure */
free_template(ATLAS_TEMPLATE * xa)1226 void free_template(ATLAS_TEMPLATE *xa)
1227 {
1228 if(xa==NULL)
1229 return;
1230 if (xa->space) free(xa->space);
1231 if (xa->template) free(xa->template);
1232 if (xa->description) free(xa->description);
1233 if (xa->comment) free(xa->comment);
1234 }
1235
1236 /* free list of atlases */
free_atlas_list(ATLAS_LIST * xal)1237 void free_atlas_list(ATLAS_LIST *xal)
1238 {
1239 int i;
1240
1241 if(xal==NULL)
1242 return;
1243 for(i=0;i<xal->natlases;i++)
1244 free_atlas(xal->atlas+i);
1245 if(xal->natlases >= 1)
1246 free(xal->atlas);
1247 free(xal);
1248 }
1249
1250 /* free an atlas dset holder */
free_adh(ATLAS_DSET_HOLDER * adh)1251 void free_adh(ATLAS_DSET_HOLDER *adh) {
1252 if (adh) free(adh);
1253 }
1254
1255 /* free an atlas structure */
free_atlas(ATLAS * xa)1256 void free_atlas(ATLAS *xa)
1257 {
1258 if(xa==NULL)
1259 return;
1260 if (xa->space) free(xa->space);
1261 if (xa->dset_name) free(xa->dset_name);
1262 if (xa->description) free(xa->description);
1263 if (xa->name) free(xa->name);
1264 if (xa->comment) free(xa->comment);
1265 if (xa->atlas_type) free(xa->atlas_type);
1266 if (xa->supp_web_info) free(xa->supp_web_info);
1267 if (xa->supp_web_type) free(xa->supp_web_type);
1268 if (xa->supp_conn_info) free(xa->supp_conn_info);
1269 if (xa->supp_conn_type) free(xa->supp_conn_type);
1270 if (xa->orient) free(xa->orient);
1271
1272 if (xa->adh) free_adh(xa->adh);
1273 }
1274
1275 /* return 1 is combined xform is identity */
is_identity_xform_list(ATLAS_XFORM_LIST * xfl,int combine)1276 int is_identity_xform_list(ATLAS_XFORM_LIST *xfl, int combine)
1277 {
1278 int i;
1279 ATLAS_XFORM_LIST *cxfl=NULL;
1280 ATLAS_XFORM *xf;
1281
1282 if(xfl==NULL) {
1283 if (wami_verb()) fprintf(stderr,"NULL transform\n");
1284 return(0);
1285 }
1286 if (combine) {
1287 if (!(cxfl = calc_xform_list(xfl))) return(0);
1288 xfl = cxfl;
1289 }
1290
1291 for(i=0;i<xfl->nxforms;i++) {
1292 xf = xfl->xform+i;
1293 if (strcmp(xf->xform_type,"Identity")) {
1294 if (cxfl) free_xform_list(cxfl);
1295 return(0);
1296 }
1297 }
1298 if (cxfl) free_xform_list(cxfl);
1299 return(1);
1300 }
1301
1302 /* return 1 if xform from src to dest is Identity */
is_identity_xform_chain(char * src,char * dest)1303 int is_identity_xform_chain(char *src, char *dest)
1304 {
1305 ATLAS_XFORM_LIST *xfl=NULL;
1306 int ans=0;
1307
1308 if (!src || !dest) return(0);
1309 if (!strcmp(src,dest)) return(1);
1310
1311 xfl = report_xform_chain(src, dest, 0);
1312 ans = is_identity_xform_list(xfl, 1);
1313 free_xform_list(xfl);
1314 return(ans);
1315 }
1316
print_xform_chain(char * src,char * dest)1317 void print_xform_chain(char *src, char *dest)
1318 {
1319 ATLAS_XFORM_LIST *xfl=report_xform_chain(src, dest,1);
1320 if (xfl) free_xform_list(xfl);
1321 return;
1322 }
1323
1324 /* print list of xforms - short form - as a chain */
print_xform_list(ATLAS_XFORM_LIST * xfl)1325 void print_xform_list(ATLAS_XFORM_LIST *xfl)
1326 {
1327 int i;
1328 ATLAS_XFORM *xf;
1329 INFO_message("----- Transform list: -------");
1330
1331 if(xfl==NULL) {
1332 fprintf(stderr,"NULL transform\n");
1333 return;
1334 }
1335 for(i=0;i<xfl->nxforms;i++) {
1336 xf = xfl->xform+i;
1337 fprintf(stderr,"%s ", xf->xform_name);
1338 if(xf->inverse)
1339 fprintf(stderr,"I");
1340 if(i==(xfl->nxforms)-1)
1341 fprintf(stderr,"\n");
1342 else
1343 fprintf(stderr," -> ");
1344 }
1345 INFO_message("");
1346 }
1347
1348 /* print list of xforms in long form */
print_all_xforms(ATLAS_XFORM_LIST * xfl)1349 void print_all_xforms(ATLAS_XFORM_LIST *xfl)
1350 {
1351 int i;
1352 ATLAS_XFORM *xf;
1353 INFO_message("----- Transform list: -------");
1354
1355 if(xfl==NULL)
1356 return;
1357 for(i=0;i<xfl->nxforms;i++) {
1358 xf = xfl->xform+i;
1359 print_xform(xf);
1360 INFO_message("-------");
1361 }
1362 INFO_message("");
1363 }
1364
1365 /* print the attributes of a transformation including its elements */
print_xform(ATLAS_XFORM * xf)1366 void print_xform(ATLAS_XFORM *xf)
1367 {
1368
1369 int i;
1370 float *xfptr;
1371
1372 fprintf(stderr, "xform: %s\n", xf->xform_name);
1373 fprintf(stderr, "xform_type: %s\n", xf->xform_type);
1374 fprintf(stderr, "xform source: %s dest: %s\n", xf->source, xf->dest);
1375 fprintf(stderr, "coord order: %s\n", xf->coord_order);
1376 fprintf(stderr, "xform dist: %f inverse: %d post: %d nelts: %d\n",
1377 xf->dist, xf->inverse, xf->post, xf->nelts);
1378 xfptr = (float *) xf->xform; /* for now use floating point values */
1379 /* transformations may require different kinds of
1380 values in the future */
1381 if(strcmp(xf->xform_type,"Affine")==0)
1382 print_affine_xform_data(xfptr);
1383 else {
1384 for (i=0;i<xf->nelts;i++)
1385 fprintf(stderr, "%f ", *xfptr++);
1386 fprintf(stderr,"\n");
1387 }
1388 }
1389
1390 /* print xform data for affine matrix in 4 columns */
print_affine_xform_data(float * xfptr)1391 void print_affine_xform_data(float *xfptr)
1392 {
1393 int i, j;
1394
1395 for (i=0;i<3;i++){
1396 for (j=0;j<4;j++)
1397 fprintf(stderr, "%f ", *xfptr++);
1398 fprintf(stderr,"\n");
1399 }
1400 fprintf(stderr,"\n");
1401 }
1402
1403 /* print list of spaces */
print_space_list(ATLAS_SPACE_LIST * xsl)1404 void print_space_list(ATLAS_SPACE_LIST *xsl)
1405 {
1406 int i;
1407 ATLAS_SPACE *xs;
1408
1409 if(xsl==NULL) {
1410 if(wami_verb() > 1)
1411 INFO_message("NULL Space list pointer, showing global list\n");
1412 xsl = get_G_space_list();
1413 }
1414 if(wami_verb() > 1)
1415 INFO_message("Space list has %d spaces\n",xsl->nspaces);
1416 INFO_message("----- List of available spaces: -------");
1417 for(i=0;i<xsl->nspaces;i++) {
1418 xs = xsl->space+i;
1419 INFO_message("%s", xs->atlas_space);
1420 }
1421 }
1422
1423 /* print list of atlases */
print_atlas_list(ATLAS_LIST * xal)1424 void print_atlas_list(ATLAS_LIST *xal)
1425 {
1426 int i;
1427 ATLAS *xa;
1428
1429 INFO_message("----- Atlas list: -------");
1430 if(xal==NULL){
1431 INFO_message("** No atlases found **");
1432 return;
1433 }
1434 for(i=0;i<xal->natlases;i++) {
1435 xa = xal->atlas+i;
1436 print_atlas(xa, 0);
1437 }
1438 }
1439
1440 /* print a single line entry description of an atlas */
print_atlas(ATLAS * xa,int level)1441 void print_atlas(ATLAS *xa, int level)
1442 {
1443 if(level==0)
1444 INFO_message("Atlas name: %s, file: %s, space: %s\n",
1445 xa->name, xa->dset_name, xa->space);
1446 else
1447 INFO_message("Atlas name: %s, file: %s, space: %s\n"
1448 "dset %p, %d sub-bricks \n"
1449 "adh %p\n",
1450 xa->name, xa->dset_name, xa->space,
1451 ATL_DSET(xa), ATL_DSET(xa) ? DSET_NVALS(ATL_DSET(xa)):-1,
1452 xa->adh);
1453 return;
1454 }
1455
1456
1457 /* print table of atlases - name, dset, description, comments */
print_atlas_table(ATLAS_LIST * xal)1458 void print_atlas_table(ATLAS_LIST *xal)
1459 {
1460 int i;
1461 ATLAS *xa;
1462
1463 INFO_message("----- Atlas list: -------");
1464 if(xal==NULL){
1465 INFO_message("** No atlases found **");
1466 return;
1467 }
1468
1469 INFO_message("Name Space Dataset Description");
1470 INFO_message("__________________________________________________________");
1471 for(i=0;i<xal->natlases;i++) {
1472 xa = xal->atlas+i;
1473 INFO_message("%-25.25s %-15.15s %s %-60s",
1474 xa->name, xa->space,
1475 xa->dset_name,
1476 ATL_DESCRIPTION_S(xa));
1477 }
1478
1479 INFO_message("\n");
1480 for(i=0;i<xal->natlases;i++) {
1481 xa = xal->atlas+i;
1482 if (ATL_COMMENT(xa)) {
1483 INFO_message("%s: %s", ATL_NAME(xa), ATL_COMMENT_S(xa));
1484 }
1485 else printf("no comment\n");
1486 }
1487 INFO_message("--------------------------");
1488 }
1489
1490 /* find dataset filename for an atlas in list of atlases */
find_atlas_dset(ATLAS_LIST * xal,char * atlasname)1491 char * find_atlas_dset(ATLAS_LIST *xal, char *atlasname)
1492 {
1493 int i;
1494 ATLAS *xa;
1495
1496 if(xal==NULL){
1497 return(NULL);
1498 }
1499
1500 // look through all the atlases
1501 for(i=0;i<xal->natlases;i++) {
1502 xa = xal->atlas+i;
1503 // is this a match?
1504 if(strcmp(xa->name, atlasname)==0){
1505 return(xa->dset_name);
1506 }
1507 }
1508 return(NULL);
1509 }
1510
1511 /* print the file name associated with an atlas */
print_atlas_dset(ATLAS_LIST * xal,char * atlasname)1512 void print_atlas_dset(ATLAS_LIST *xal, char *atlasname)
1513 {
1514 char *dsetname;
1515
1516 dsetname = find_atlas_dset(xal, atlasname);
1517 if(dsetname) printf("%s\n", dsetname);
1518 }
1519
1520
1521 /* print the comment for an atlas - may span multiple lines with '\n' */
print_atlas_comment(ATLAS * xa)1522 void print_atlas_comment(ATLAS *xa)
1523 {
1524 if((xa) && ATL_COMMENT(xa))
1525 INFO_message("%s", ATL_COMMENT(xa));
1526 }
1527
1528 /* print the atlas_type for an atlas */
print_atlas_type(ATLAS * xa)1529 void print_atlas_type(ATLAS *xa)
1530 {
1531 if((xa) && ATL_TYPE(xa))
1532 INFO_message("%s", ATL_TYPE(xa));
1533 }
1534
1535 /* print the atlas supplemental information for an atlas */
print_atlas_supp_web_info(ATLAS * xa)1536 void print_atlas_supp_web_info(ATLAS *xa)
1537 {
1538 if((xa) && ATL_SUPP_WEB_INFO(xa)) {
1539 INFO_message("%sroiname%s", ATL_SUPP_WEB_INFO_S(xa),
1540 ATL_SUPP_WEB_TYPE_S(xa));
1541 }
1542 if((xa) && ATL_SUPP_CONN_INFO(xa)) {
1543 INFO_message("%sroiname%s", ATL_SUPP_CONN_INFO_S(xa),
1544 ATL_SUPP_CONN_TYPE_S(xa));
1545 }
1546
1547 }
1548
1549 /* print the coordinate orientation for an atlas */
print_atlas_orient(ATLAS * xa)1550 void print_atlas_orient(ATLAS *xa)
1551 {
1552 if((xa) && ATL_ORIENT(xa))
1553 INFO_message("%s", ATL_ORIENT(xa));
1554 }
1555
1556 /* print table of atlases - name, dset, description, comments */
print_point_lists(ATLAS_LIST * xal)1557 void print_point_lists(ATLAS_LIST *xal)
1558 {
1559 int i;
1560 ATLAS *xa;
1561 ATLAS_POINT_LIST *apl;
1562
1563 INFO_message("----- Atlas point lists: -------");
1564 if(xal==NULL){
1565 INFO_message("** No atlases found **");
1566 return;
1567 }
1568
1569 for(i=0;i<xal->natlases;i++) {
1570 xa = xal->atlas+i;
1571 INFO_message("Atlas name : %-25.25s, Dataset: %-54.54s",
1572 xa->name,xa->dset_name);
1573 apl = atlas_point_list(xa->name);
1574 if(apl)
1575 print_atlas_point_list(apl);
1576 else{
1577 if(ATL_WEB_TYPE(xa))
1578 INFO_message("web-based atlas. No local point list");
1579 else
1580 INFO_message("**** No point list. Atlas needs repair!");
1581 }
1582 INFO_message(
1583 "__________________________________________________________");
1584 }
1585
1586 INFO_message("\n");
1587 for(i=0;i<xal->natlases;i++) {
1588 xa = xal->atlas+i;
1589 if (ATL_COMMENT(xa)) {
1590 INFO_message("%s: %s", ATL_NAME(xa), ATL_COMMENT_S(xa));
1591 }
1592 }
1593 INFO_message("--------------------------");
1594 }
1595
1596 /* return max length of a label in an atlas point list */
1597 int
atlas_max_label_length(ATLAS_POINT * ap,int n_points)1598 atlas_max_label_length(ATLAS_POINT *ap, int n_points)
1599 {
1600 int i,len,maxlen=0;
1601
1602 if(ap==NULL)
1603 return(0);
1604
1605 for(i=0;i<n_points;i++) {
1606 /* name (+long name?) */
1607 len = strlen(Atlas_name_choice(&ap[i]));
1608 if(len>maxlen) maxlen = len;
1609 }
1610
1611 return(maxlen);
1612 }
1613
1614 /* mark if any atlas points have a non-zero level associated with them */
1615 int
atlas_level(ATLAS_POINT * ap,int n_points)1616 atlas_level(ATLAS_POINT *ap, int n_points)
1617 {
1618 int i;
1619
1620 if(ap==NULL)
1621 return(0);
1622
1623 for(i=0;i<n_points;i++) {
1624 if(ap[i].tdlev) return(1);
1625 }
1626
1627 return(0);
1628 }
1629
1630 /* print info about an atlas coordinate */
print_atlas_coord(ATLAS_COORD ac)1631 void print_atlas_coord (ATLAS_COORD ac)
1632 {
1633 INFO_message("----- Atlas Coord: -------");
1634 INFO_message("%f %f %f (%s), space_name %s%s\n",
1635 ac.x, ac.y, ac.z, ac.orcode, ac.space_name,
1636 is_known_coord_space(ac.space_name) ? "":"Lost in space");
1637 }
1638
is_known_coord_space(char * space_name)1639 int is_known_coord_space(char *space_name) {
1640 if (!space_name || space_name[0]=='\0' ||
1641 !strcmp(space_name,"Unknown")) return(0);
1642 if (find_atlas_space_index(space_name) >= 0) return(1);
1643 return(0);
1644 }
1645
1646 /* print list of templates */
print_template_list(ATLAS_TEMPLATE_LIST * xtl)1647 void print_template_list(ATLAS_TEMPLATE_LIST *xtl)
1648 {
1649 int i, templen;
1650 ATLAS_TEMPLATE *xt;
1651 char *tempstr, *tempstr2=NULL;
1652 INFO_message("----- Template list: -------");
1653 if(xtl==NULL)
1654 return;
1655 for(i=0;i<xtl->ntemplates;i++) {
1656 xt = xtl->atlas_template+i;
1657 if(xt->description) {
1658 tempstr2 = ATL_DESCRIPTION_S(xt);
1659 templen = strlen(xt->template)+strlen(tempstr2) + 3;
1660 tempstr = (char *)calloc( templen, sizeof(char));
1661 sprintf(tempstr, "%s: %s", xt->template, ATL_DESCRIPTION_S(xt));
1662 show_wrapping_line(tempstr,"", 0, stdout);
1663 free(tempstr);
1664 }
1665 else {
1666 show_wrapping_line(xt->template,"", 0, stdout);
1667 }
1668 if(xt->comment){
1669 /* show_wrapping_line("Comment:\n","", 0, stdout);*/
1670 show_wrapping_line(ATL_COMMENT(xt)," ", 0, stdout);
1671 }
1672 }
1673 }
1674
1675
1676 /* print the comment for a template - may span multiple lines with '\n' */
print_template_comment(ATLAS_TEMPLATE * xa)1677 void print_template_comment(ATLAS_TEMPLATE *xa)
1678 {
1679 if((xa) && ATL_COMMENT(xa))
1680 INFO_message("%s", ATL_COMMENT(xa));
1681 }
1682
1683
1684 /* apply line indent per line, if we exceed MAX_LINE_CHARS, wrap */
1685 /* from afni_history.c - rickr, moved by drg */
show_wrapping_line(char * str,char * prefix,int indent,FILE * fp)1686 int show_wrapping_line(char * str, char * prefix, int indent, FILE * fp)
1687 {
1688 int c, cline, len;
1689
1690 if( !str ) return 0;
1691
1692 if( prefix ) fputs(prefix, fp);
1693
1694 len = strlen(str);
1695 if( len < 2 ) return 1;
1696
1697 if( str[len-1] == '\n' ) len--; /* ignore trailing newline */
1698
1699 cline = 0;
1700 for( c = 0; c < len; c++ ) {
1701 if( str[c] == '\n' ) { /* print newline and indent */
1702 fputc('\n', fp);
1703 fprintf(fp, "%*s", indent, "");
1704 cline = 0;
1705 continue;
1706 } else if ( cline > MAX_WAMI_LINE_CHARS ) { /* fix, and continue */
1707 fputc('\n', fp);
1708 fprintf(fp, "%*s", indent, "");
1709 cline = 0;
1710 }
1711 fputc(str[c], fp);
1712 cline++;
1713 }
1714
1715 fprintf(fp,"\n");
1716
1717 return 0;
1718 }
1719
1720
1721 /* calculate list of xforms */
calc_xform_list(ATLAS_XFORM_LIST * xfl)1722 ATLAS_XFORM_LIST *calc_xform_list(ATLAS_XFORM_LIST *xfl)
1723 {
1724 int i, nxf, sl1, sl2, cc;
1725 ATLAS_XFORM *xf, *xf2, *xf3, *oldxfptr=NULL;
1726 char *source, *dest;
1727 ATLAS_XFORM_LIST *cxfl;
1728
1729 if(wami_verb() > 1)
1730 printf("calculating xform list\n");
1731 if(xfl==NULL)
1732 return(NULL);
1733 nxf = (xfl->nxforms) - 1;
1734
1735 /* make condensed transformation list structure */
1736 cxfl = (ATLAS_XFORM_LIST *)calloc(1,sizeof(ATLAS_XFORM_LIST));
1737 if(cxfl==NULL)
1738 ERROR_exit("Could not allocate space for condensed xform list\n");
1739 cxfl->xform = (ATLAS_XFORM *) calloc((xfl->nxforms),sizeof(ATLAS_XFORM));
1740 if(cxfl->xform==NULL)
1741 ERROR_exit("Could not allocate space for condensed xform list xforms\n");
1742
1743 cxfl->nxforms = 0;
1744 if(wami_verb() > 1)
1745 printf("starting to combine xforms\n");
1746 /* only one xform, just go home */
1747 if(xfl->nxforms==1){
1748 if(wami_verb() > 1)
1749 printf("only 1 xform\n");
1750 cxfl->nxforms = 1;
1751 cc = copy_xform(xfl->xform, cxfl->xform);
1752 if(cc!=0) {
1753 ERROR_exit("Could not copy only xform for condensed xform list");
1754 }
1755 if(cxfl->xform->inverse == 1){ /* check for inverse, even if only one */
1756 xf = cxfl->xform;
1757 invert_xform(cxfl->xform);
1758 source = nifti_strdup(xf->dest);
1759 dest = nifti_strdup(xf->source);
1760 free(xf->xform_name);
1761 free(xf->source); free(xf->dest);
1762 xf->source = source;
1763 xf->dest = dest;
1764 sl1 = strlen(xf->source); sl2 = strlen(xf->dest);
1765 xf->xform_name = (char *) calloc((sl1+sl2+3),sizeof(char));
1766 sprintf(xf->xform_name, "%s::%s", xf->source, xf->dest);
1767 }
1768 return(cxfl);
1769 }
1770 /* do calculations on xforms a pair at a time */
1771 xf = xfl->xform;
1772 for(i=0;i<nxf;i++) {
1773 if(wami_verb() > 1)
1774 printf("xf %d with xf %d\n", i, i+1);
1775
1776 xf2 = xfl->xform+i+1;
1777
1778 if(xf2->inverse)
1779 dest = nifti_strdup(xf2->source);
1780 else
1781 dest = nifti_strdup(xf2->dest);
1782
1783 if(xf->inverse)
1784 source = nifti_strdup(xf->dest);
1785 else
1786 source = nifti_strdup(xf->source);
1787
1788 if(wami_verb() > 1)
1789 INFO_message("Multiplying %s type with %s type in chain\n", xf->xform_type,
1790 xf2->xform_type);
1791
1792 xf3 = calc_xf(xf,xf2); /* calculate xforms a pair of time */
1793 if(xf3) {
1794 free(xf3->xform_name);
1795 free(xf3->source); free(xf3->dest);
1796
1797 xf3->source = source;
1798 xf3->dest = dest;
1799 sl1 = strlen(xf3->source); sl2 = strlen(xf3->dest);
1800 xf3->xform_name = (char *) calloc((sl1+sl2+3),sizeof(char));
1801 sprintf(xf3->xform_name, "%s::%s", xf3->source, xf3->dest);
1802
1803 if(i==(nxf-1)){
1804 if(wami_verb() > 1)
1805 printf(
1806 "On last xform, copying last combined xform"
1807 " to combined xform list\n");
1808 cc = copy_xform(xf3, cxfl->xform+(cxfl->nxforms));
1809 cxfl->nxforms++;
1810 if(wami_verb() > 1){
1811 print_xform(xf3);
1812 xf = xf3;
1813 print_xform(xf);
1814 }
1815 }
1816 else{
1817 if(wami_verb() > 1)
1818 printf("could combine xform %d with %d\n", i, i+1);
1819 xf = xf3; cc = 0;
1820 /* cc = copy_xform(xf3, xf); */ /* use combined xform for next pair */
1821 if(wami_verb() > 1)
1822 print_xform(xf);
1823
1824 }
1825
1826 }
1827 else {
1828 if(wami_verb() > 1)
1829 printf("could not calculate this combination of xforms "
1830 "- adding to chain\n");
1831 cc = copy_xform(xf, cxfl->xform+(cxfl->nxforms));
1832 cxfl->nxforms++;
1833 /* if on last combination, add the last one too */
1834 if(i==nxf-1){
1835 cc = copy_xform(xf2, cxfl->xform+(cxfl->nxforms));
1836 cxfl->nxforms++;
1837 }
1838
1839 if((cc==0)&&(i<nxf-1))
1840 xf = xf2; cc = 0;
1841 /*copy_xform(xf2, xf); */ /* update start xform for next pair */
1842 }
1843
1844 if(i>0) /* free the temporary xform intermediate */
1845 free_xform(oldxfptr);
1846 oldxfptr = xf3;
1847
1848 if(cc!=0) {
1849 ERROR_exit("Could not copy a xform for condensed xform list");
1850 }
1851
1852 }
1853
1854
1855 return(cxfl);
1856 }
1857
1858 /* calculate product for pair of xforms */
1859 /* return xform product - allocating space even if copy
1860 If xf and xf2 are both affine transformations,
1861 return xf2 * xf */
1862 ATLAS_XFORM *
calc_xf(ATLAS_XFORM * xf,ATLAS_XFORM * xf2)1863 calc_xf(ATLAS_XFORM *xf, ATLAS_XFORM *xf2)
1864 {
1865 ATLAS_XFORM *xf3;
1866 int cc;
1867
1868 /* if(wami_verb() > 1)
1869 INFO_message("Multiplying %s type with %s type\n", xf->xform_type,
1870 xf2->xform_type);
1871 */
1872
1873 xf3 = calloc(1,sizeof(ATLAS_XFORM));
1874 if(xf3==NULL)
1875 return(NULL);
1876 invert_xform(xf); /* possibly need to invert transform */
1877 invert_xform(xf2);
1878
1879
1880 /* check for identity transformations - simplest */
1881 if(strcmp(xf->xform_type,"Identity")==0){
1882 cc = copy_xform(xf2,xf3);
1883 if(cc!=0) {
1884 return(NULL);
1885 }
1886 else return(xf3);
1887 }
1888
1889 if(strcmp(xf2->xform_type,"Identity")==0){
1890 cc = copy_xform(xf,xf3);
1891 if(cc!=0) {
1892 return(NULL);
1893 }
1894 else return(xf3);
1895 }
1896
1897
1898 if(wami_verb() > 1)
1899 INFO_message("Multiplying %s type with %s type\n", xf->xform_type,
1900 xf2->xform_type);
1901 if(strcmp(xf->xform_type,"Affine")==0){
1902 if(strcmp(xf2->xform_type,"Affine")==0){
1903 cc = affine_mult(xf,xf2,xf3);
1904 if(wami_verb() > 1)
1905 INFO_message("combining two affine matrices\n");
1906 if(cc!=0) {
1907 if(wami_verb() > 1)
1908 INFO_message("could not combine two affine matrices\n");
1909
1910 return(NULL);
1911 }
1912 else return(xf3);
1913 }
1914 if(strcmp(xf2->xform_type,"2-piece")==0){
1915 cc = affine_2piece_mult(xf,xf2,xf3,0);
1916 if(cc!=0) {
1917 return(NULL);
1918 }
1919 else return(xf3);
1920 }
1921 if(strcmp(xf2->xform_type,"12-piece")==0){
1922 cc = affine_12piece_mult(xf,xf2,xf3,0);
1923 if(cc!=0) {
1924 return(NULL);
1925 }
1926 else return(xf3);
1927 }
1928 }
1929
1930 if(strcmp(xf->xform_type,"2-piece")==0){
1931 if(strcmp(xf2->xform_type,"Affine")==0){
1932 cc = affine_2piece_mult(xf,xf2,xf3,-1);
1933 if(cc!=0) {
1934 return(NULL);
1935 }
1936 else return(xf3);
1937 }
1938 if(strcmp(xf2->xform_type,"2-piece")==0){
1939 cc = x2piece_2piece_mult(xf,xf2,xf3);
1940 if(cc!=0) {
1941 return(NULL);
1942 }
1943 else return(xf3);
1944 }
1945 if(strcmp(xf2->xform_type,"12-piece")==0){
1946 cc = x2piece_12piece_mult(xf,xf2,xf3,0);
1947 if(cc!=0) {
1948 return(NULL);
1949 }
1950 else return(xf3);
1951 }
1952 }
1953
1954 if(strcmp(xf->xform_type,"12-piece")==0){
1955 if(strcmp(xf2->xform_type,"Affine")==0){
1956 cc = affine_12piece_mult(xf,xf2,xf3,-1);
1957 if(cc!=0) {
1958 return(NULL);
1959 }
1960 else return(xf3);
1961 }
1962 if(strcmp(xf2->xform_type,"2-piece")==0){
1963 cc = x2piece_12piece_mult(xf,xf2,xf3,-1);
1964 if(cc!=0) {
1965 return(NULL);
1966 }
1967 else return(xf3);
1968 }
1969 if(strcmp(xf2->xform_type,"12-piece")==0){
1970 cc = x12piece_12piece_mult(xf,xf2,xf3);
1971 if(cc!=0) {
1972 return(NULL);
1973 }
1974 else return(xf3);
1975 }
1976 }
1977
1978 if(wami_verb())
1979 INFO_message("AFNI doesn't know how to combine these transforms\n"
1980 "Using the transformations sequentially");
1981
1982 return(NULL);
1983 }
1984
1985 /* invert transformation - general form */
1986 int
invert_xform(ATLAS_XFORM * xf)1987 invert_xform(ATLAS_XFORM *xf)
1988 {
1989 int cc =1;
1990
1991 if(xf->inverse==0) return(0);
1992
1993 xf->inverse = 0;
1994
1995 if(strcmp(xf->xform_type,"Identity")==0){
1996 return(0); /* inverse of identity is identity */
1997 }
1998
1999 if(strcmp(xf->xform_type,"Affine")==0){
2000 cc = invert_affine(xf);
2001 }
2002
2003 /* for 12-piece, this just flips inverse setting 0->1*/
2004 if(strcmp(xf->xform_type,"12-piece")==0){
2005 cc = invert_12piece(xf);
2006 }
2007 if(strcmp(xf->xform_type,"2-piece")==0){
2008 cc = invert_2piece(xf);
2009 }
2010 if(strcmp(xf->xform_type,"brett_mni2tt")==0){
2011 cc = invert_brett(xf);
2012 }
2013
2014
2015 return(cc);
2016 }
2017
2018 /* invert an affine matrix - do in place */
2019 /* return error code - 0 for no error */
2020 int
invert_affine(ATLAS_XFORM * xf)2021 invert_affine(ATLAS_XFORM *xf)
2022 {
2023 int i, j;
2024 matrix tempmat, invmat;
2025 float *xfptr;
2026 ENTRY("invert_affine");
2027
2028 if ( !xf || !xf->xform) RETURN(1);
2029
2030 matrix_initialize (&tempmat);
2031 matrix_create(4,4,&tempmat);
2032 xfptr = (float *) xf->xform;
2033 for(i=0;i<3;i++)
2034 for(j=0;j<4;j++)
2035 tempmat.elts[i][j] = (double) *xfptr++; /* recast float to double*/
2036 tempmat.elts[3][0] = tempmat.elts[3][1] = tempmat.elts[3][2] = 0.0;
2037 tempmat.elts[3][3] = 1.0;
2038 matrix_initialize (&invmat);
2039 matrix_inverse(tempmat, &invmat);
2040
2041 xfptr = (float *) xf->xform;
2042 for(i=0;i<3;i++)
2043 for(j=0;j<4;j++)
2044 *xfptr++ = (float) invmat.elts[i][j];
2045
2046 matrix_destroy(&invmat);
2047 matrix_destroy(&tempmat);
2048
2049 RETURN(0);
2050 }
2051
2052 /* invert a 12 piece matrix - do in place */
2053 /* don't actually invert, just mark for backward vector use */
2054 /* return error code - 0 for no error */
2055 int
invert_12piece(ATLAS_XFORM * xf)2056 invert_12piece(ATLAS_XFORM *xf)
2057 {
2058 xf->inverse = !(xf->inverse);
2059
2060 return(0);
2061 }
2062
2063
2064 /* invert a 2 piece matrix - do in place */
2065 /* return error code - 0 for no error */
2066 int
invert_2piece(ATLAS_XFORM * xf)2067 invert_2piece(ATLAS_XFORM *xf)
2068 {
2069 /* int cc;
2070 */
2071 return(1);
2072 }
2073
2074 /* invert a brett transform - do in place */
2075 /* return error code - 0 for no error */
2076 int
invert_brett(ATLAS_XFORM * xf)2077 invert_brett(ATLAS_XFORM *xf)
2078 {
2079 xf->inverse = !(xf->inverse);
2080
2081 return(0);
2082 }
2083
2084 /* multiply affine transformations */
2085 int
affine_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3)2086 affine_mult(ATLAS_XFORM *xf, ATLAS_XFORM *xf2, ATLAS_XFORM *xf3)
2087 {
2088 int cc, i, j;
2089 matrix sm1, sm2, sm3, sm4;
2090 float *xfptr, *xfptr2;
2091 cc = copy_xform(xf,xf3); /* allocate xform structure */
2092 if (cc!=0)
2093 return(1);
2094
2095 matrix_initialize(&sm1);
2096 matrix_initialize(&sm2);
2097 matrix_initialize(&sm3);
2098 matrix_create(4,4,&sm1);
2099 matrix_create(4,4,&sm2);
2100
2101 xfptr = (float *) xf->xform;
2102 xfptr2 = (float *) xf2->xform;
2103 for(i=0;i<3;i++)
2104 for(j=0;j<4;j++) {
2105 sm1.elts[i][j] = (double) *xfptr++;
2106 sm2.elts[i][j] = (double) *xfptr2++;
2107 }
2108
2109 /* fill in last line with 0's and 1's to make invertible 4x4 */
2110 sm1.elts[3][0] = sm1.elts[3][1] = sm1.elts[3][2] = 0.0;
2111 sm1.elts[3][3] = 1.0;
2112 sm2.elts[3][0] = sm2.elts[3][1] = sm2.elts[3][2] = 0.0;
2113 sm2.elts[3][3] = 1.0;
2114
2115 /* if xforms are different coordinate order (RAI/LPI)
2116 apply -1,-1,1 diagonal matrix to first matrix */
2117 /* may want to expand for all other coord orders */
2118 if(strcmp(xf->coord_order, xf2->coord_order)!=0) {
2119 matrix_initialize(&sm4);
2120 matrix_identity(4,&sm4);
2121 sm4.elts[0][0] = -1.0;
2122 sm4.elts[1][1] = -1.0;
2123 matrix_multiply(sm1, sm4, &sm3);
2124 matrix_multiply(sm4, sm3, &sm1);
2125 matrix_destroy(&sm4);
2126 if(xf3->coord_order) free(xf3->coord_order);
2127 xf3->coord_order = nifti_strdup(xf2->coord_order);
2128 }
2129
2130 /* multiply destination direction first, then source s3->s2 x s2->s1 */
2131 /* mod - drg 02/27/20011 - was backwards, sm1*sm2 - yikes! */
2132 matrix_multiply(sm2, sm1, &sm3);
2133
2134 xfptr = (float *) xf3->xform;
2135 for(i=0;i<3;i++)
2136 for(j=0;j<4;j++) {
2137 *xfptr++ = (float) sm3.elts[i][j];
2138 }
2139
2140 matrix_destroy(&sm1);
2141 matrix_destroy(&sm2);
2142 matrix_destroy(&sm3);
2143
2144 return(0);
2145 }
2146
2147 /* multiply affine transformation by 2-piece transform */
2148 int
affine_2piece_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3,int dir)2149 affine_2piece_mult(ATLAS_XFORM *xf, ATLAS_XFORM *xf2, ATLAS_XFORM *xf3, int dir)
2150 {
2151 int cc;
2152
2153 return(1); /* can't do this yet */
2154 if(dir)
2155 cc = copy_xform(xf2,xf3);
2156 else
2157 cc = copy_xform(xf,xf3); /* allocate xform structure */
2158 if (cc!=0)
2159 return(1);
2160 return(0);
2161 }
2162
2163
2164 /* multiply affine transformation by 12-piece transform */
2165 int
affine_12piece_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3,int dir)2166 affine_12piece_mult(ATLAS_XFORM *xf,
2167 ATLAS_XFORM *xf2, ATLAS_XFORM *xf3, int dir)
2168 {
2169 int cc;
2170
2171 return(1);
2172
2173 if(dir)
2174 cc = copy_xform(xf2,xf3);
2175 else
2176 cc = copy_xform(xf,xf3); /* allocate xform structure */
2177 if (cc!=0)
2178 return(1);
2179 return(0);
2180 }
2181
2182 /* multiply two 2-piece affine transformations */
2183 int
x2piece_2piece_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3)2184 x2piece_2piece_mult(ATLAS_XFORM *xf,
2185 ATLAS_XFORM *xf2, ATLAS_XFORM *xf3)
2186 {
2187 int cc;
2188
2189 return(1);
2190 cc = copy_xform(xf,xf3); /* allocate xform structure */
2191 if (cc!=0)
2192 return(1);
2193 return(0);
2194 }
2195
2196 /* multiply a 2-piece and a 12-piece affine transformation */
2197 int
x2piece_12piece_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3,int dir)2198 x2piece_12piece_mult(ATLAS_XFORM *xf,
2199 ATLAS_XFORM *xf2, ATLAS_XFORM *xf3, int dir)
2200 {
2201 int cc;
2202
2203 return(1);
2204 if(dir)
2205 cc = copy_xform(xf2,xf3);
2206 else
2207 cc = copy_xform(xf,xf3); /* allocate xform structure */
2208 if (cc!=0)
2209 return(1);
2210 return(0);
2211 }
2212
2213
2214 /* multiply two 12-piece affine transformations */
2215 int
x12piece_12piece_mult(ATLAS_XFORM * xf,ATLAS_XFORM * xf2,ATLAS_XFORM * xf3)2216 x12piece_12piece_mult(ATLAS_XFORM *xf,
2217 ATLAS_XFORM *xf2, ATLAS_XFORM *xf3)
2218 {
2219 int cc;
2220
2221 return(1);
2222 cc = copy_xform(xf,xf3); /* allocate xform structure */
2223 if (cc!=0)
2224 return(1);
2225 return(0);
2226 }
2227
2228
2229 /* apply xform to xyz */
2230 int
apply_xform_general(ATLAS_XFORM * xf,float x,float y,float z,float * xout,float * yout,float * zout)2231 apply_xform_general(ATLAS_XFORM *xf, float x,float y,float z,
2232 float *xout, float *yout, float *zout)
2233 {
2234 int xgc = 1;
2235
2236 invert_xform(xf); /* possibly need to invert transform */
2237
2238 if(strcmp(xf->xform_type,"Affine")==0){
2239 xgc = apply_xform_affine(xf, x, y, z, xout, yout, zout);
2240 }
2241 if(strcmp(xf->xform_type,"2-piece")==0){
2242 xgc = apply_xform_2piece(xf, x, y, z, xout, yout, zout);
2243 }
2244
2245 if(strcmp(xf->xform_type,"brett_tt2mni")==0){
2246 if(!xf->inverse)
2247 xgc = apply_xform_brett_tt2mni(x, y, z, xout, yout, zout);
2248 else
2249 xgc = apply_xform_brett_mni2tt(x, y, z, xout, yout, zout);
2250 }
2251
2252 if(strcmp(xf->xform_type,"brett_mni2tt")==0){
2253 if(xf->inverse)
2254 xgc = apply_xform_brett_tt2mni(x, y, z, xout, yout, zout);
2255 else
2256 xgc = apply_xform_brett_mni2tt(x, y, z, xout, yout, zout);
2257 }
2258
2259 if(strcmp(xf->xform_type,"12-piece")==0){
2260 xgc = apply_xform_12piece(xf, x, y, z, xout, yout, zout);
2261 }
2262 if(strcmp(xf->xform_type,"Identity")==0){
2263 *xout = x; *yout = y; *zout = z; xgc = 0;
2264 }
2265
2266 if(wami_verb() > 2)
2267 INFO_message("xform RAI out x y z %f %f %f", *xout,*yout,*zout);
2268
2269 return(xgc);
2270 }
2271
2272
2273 /* apply xform chain to xyz */
2274 int
apply_xform_chain(ATLAS_XFORM_LIST * xfl,float x,float y,float z,float * xout,float * yout,float * zout)2275 apply_xform_chain(ATLAS_XFORM_LIST *xfl, float x, float y, float z,
2276 float *xout, float *yout, float *zout)
2277 {
2278
2279 int i, nxf, xgc;
2280 float xxout=0.0, yyout=0.0, zzout=0.0;
2281 ATLAS_XFORM * xf;
2282
2283 *xout = xxout;
2284 *yout = yyout;
2285 *zout = zzout;
2286 if(!xfl || !xfl->xform) {
2287 return(0);
2288 }
2289 nxf = xfl->nxforms;
2290 if(nxf==0) return(0);
2291
2292 for(i=0;i<nxf;i++) {
2293 xf = xfl->xform+i;
2294 xgc = apply_xform_general(xf, x, y, z, &xxout, &yyout, &zzout);
2295 if(xgc==0) {
2296 x = xxout;
2297 y = yyout;
2298 z = zzout;
2299 }
2300 else {
2301 WARNING_message("Could not transform between spaces");
2302 return(-1);
2303 }
2304 }
2305
2306 *xout = xxout;
2307 *yout = yyout;
2308 *zout = zzout;
2309
2310 return(0);
2311 }
2312
2313 /* apply the forward affine transformation to the xyz */
2314 int
apply_xform_affine(ATLAS_XFORM * xf,float x,float y,float z,float * xout,float * yout,float * zout)2315 apply_xform_affine(ATLAS_XFORM *xf, float x, float y, float z, \
2316 float *xout, float *yout, float *zout)
2317 {
2318 float *xfptr;
2319
2320 if(xf->xform==NULL) return(1);
2321
2322 /* should expand coord_order to other orientations */
2323 if(strcmp(xf->coord_order,"lpi") == 0){
2324 x = -x; y =-y;
2325 }
2326
2327
2328 xfptr = (float *) xf->xform;
2329
2330 xfptr = (float *) xf->xform;
2331
2332 *xout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2333 xfptr += 4;
2334 *yout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2335 xfptr += 4;
2336 *zout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2337
2338 if(strcmp(xf->coord_order,"lpi") == 0){
2339 *xout = -(*xout); *yout= -(*yout);
2340 }
2341
2342
2343 return(0);
2344 }
2345
2346 /* apply the forward 2-piece transformation to the xyz coordinate
2347 2-piece transformations apply one of two affine transformations
2348 depending on some dividing plane defined by x,y or z equal to
2349 a particular value. This dividing plane can be defined on the data
2350 before the transformation or it may be defined 'post' transformation.
2351 If a post transformation, the second xform applies an alternate
2352 xform to the coordinate transformed by the first affine xform.
2353 This diagram shows the various possibilities:
2354
2355 pre-xform (post=0)
2356 any x,y,z < lx,ly,lz, apply 2nd xform
2357 else apply 1st xform
2358 post-xform (post=1)
2359 apply 1st xform to get x',y',z'
2360 any x',y',z' < lx,ly,lz, apply 2nd xform to x',y',z'
2361
2362 */
2363 int
apply_xform_2piece(ATLAS_XFORM * xf,float x,float y,float z,float * xout,float * yout,float * zout)2364 apply_xform_2piece(ATLAS_XFORM *xf, float x, float y, float z,
2365 float *xout, float *yout, float *zout)
2366 {
2367 float *xfptr;
2368 float lx,ly,lz;
2369 int apply_2nd;
2370
2371 /* brett transform - tta to mni
2372 1.0101 0 0 0
2373 0 1.02962 -0.05154 0
2374 0 0.05434 1.08554 0
2375
2376 1 0 0 0
2377 0 1 0 0
2378 0 0 1.09523 0
2379 */
2380 if(!xf || xf->xform==NULL) return(1);
2381
2382 xfptr = xf->xform;
2383 /* change input coords to lpi if xform defined that way (RAI input)*/
2384 if(strcmp(xf->coord_order,"lpi") == 0){
2385 x = -x; y =-y;
2386 }
2387
2388 /* x,y,z limits are the first three elements of the xform data */
2389 lx = *xfptr++; ly = *xfptr++; lz = *xfptr++;
2390 /* If this is a pre transformation (post=0), check the x,y,z limits.
2391 If any xyz > limit, use the second transformation matrix */
2392 if(!xf->post){
2393 apply_2nd = 0;
2394 if(lx > -9998) {
2395 if(x<lx) apply_2nd = 1;
2396 }
2397 if(ly > -9998) {
2398 if(y<ly) apply_2nd = 1;
2399 }
2400 if(lz > -9998) {
2401 if(z<lz) apply_2nd = 1;
2402 }
2403 if(apply_2nd)
2404 xfptr += 12;
2405 }
2406
2407 *xout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2408 xfptr += 4;
2409 *yout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2410 xfptr += 4;
2411 *zout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2412
2413 if(xf->post){
2414 apply_2nd = 0;
2415 if(lx > -9998) {
2416 if(x<lx) apply_2nd = 1;
2417 }
2418 if(ly > -9998) {
2419 if(y<ly) apply_2nd = 1;
2420 }
2421 if(lz > -9998) {
2422 if(z<lz) apply_2nd = 1;
2423 }
2424 if(apply_2nd) {
2425 x = *xout; y = *yout; z = *zout; /* if not using a concatenated matrix*/
2426 xfptr += 4;
2427 *xout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2428 xfptr += 4;
2429 *yout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2430 xfptr += 4;
2431 *zout = (*xfptr * x) + (*(xfptr+1) * y) + (*(xfptr+2) * z) + *(xfptr+3);
2432 }
2433 }
2434 if(strcmp(xf->coord_order,"lpi") == 0){
2435 *xout = -(*xout); *yout= -(*yout);
2436 }
2437
2438 return(0);
2439 }
2440
2441 /* apply the forward 12-piece transformation to the xyz */
2442 int
apply_xform_12piece(ATLAS_XFORM * xf,float x,float y,float z,float * xout,float * yout,float * zout)2443 apply_xform_12piece(ATLAS_XFORM *xf, float x, float y, float z, \
2444 float *xout, float *yout, float *zout)
2445 {
2446 THD_talairach_12_warp *ww;
2447 int iw, ioff;
2448 THD_fvec3 tv, mv;
2449 float tx, ty, tz;
2450 char *wptr;
2451
2452 if(xf->xform==NULL) return(1);
2453 LOAD_FVEC3( mv , x,y,z ) ;
2454
2455 /* load the transform */
2456 ww = myRwcNew( THD_talairach_12_warp ) ;
2457 ww->type = WARP_TALAIRACH_12_TYPE;
2458 ww->resam_type = 0;
2459
2460 for (iw=0; iw < 12; ++iw) {
2461 /* 12 piece tranformations stored in 12x30 float blocks */
2462 /* each 30 float block is first 3x3 forward, 3x3 backward,
2463 3 vector-forward, 3 vector-backward,
2464 bottom xyz corner, top xyz corner */
2465 #if 0
2466 THD_fvec3 *fvptr;
2467 fptr = (float *) xf->xform+iw*MAPPING_LINEAR_FSIZE;
2468 memcpy(&ww->warp[iw].mfor, fptr, 9*sizeof(float)); fptr += 9;
2469 memcpy(&ww->warp[iw].mbac, fptr, 9*sizeof(float)); fptr += 9;
2470 memcpy(&ww->warp[iw].bvec, fptr, 3*sizeof(float)); fptr += 3;
2471 memcpy(&ww->warp[iw].svec, fptr, 3*sizeof(float)); fptr += 3;
2472 memcpy(&ww->warp[iw].bot, fptr, 3*sizeof(float)); fptr += 3;
2473 memcpy(&ww->warp[iw].top, fptr, 3*sizeof(float));
2474 #endif
2475
2476 ww->warp[iw].type = MAPPING_LINEAR_TYPE ;
2477
2478 ioff = iw * MAPPING_LINEAR_FSIZE ; /* increments by 12=4x3 */
2479 wptr = (char *) xf->xform + iw*MAPPING_LINEAR_FSIZE*sizeof(float);
2480 COPY_INTO_STRUCT( ww->warp[iw] ,
2481 MAPPING_LINEAR_FSTART ,
2482 float ,
2483 wptr,
2484 MAPPING_LINEAR_FSIZE ) ;
2485
2486
2487 }
2488
2489 if (!ww) {
2490 ERROR_message("Failed to form built-in warp.");
2491 return(1);
2492 } else {
2493 /*THD_fvec3 AFNI_forward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )*/
2494 if (! xf->inverse )
2495 tv = AFNI_forward_warp_vector((THD_warp *)ww, mv);
2496 else
2497 tv = AFNI_backward_warp_vector((THD_warp *)ww, mv);
2498 }
2499
2500 UNLOAD_FVEC3(tv, tx, ty, tz);
2501 *xout = tx; *yout = ty; *zout =tz;
2502
2503 free(ww);
2504
2505 return(0);
2506 }
2507
2508 /* apply Brett transform to transform from TT to MNI (Talairach to MNI)*/
2509 /* special case of 2-piece transform and maybe the only one we will ever use */
2510 int
apply_xform_brett_tt2mni(float x,float y,float z,float * xout,float * yout,float * zout)2511 apply_xform_brett_tt2mni(float x, float y, float z, \
2512 float *xout, float *yout, float *zout)
2513 {
2514
2515 /* brett transform - tta to mni
2516 1.0101 0 0 0
2517 0 1.02962 -0.05154 0
2518 0 0.05434 1.08554 0
2519
2520 1.0101 0 0 0
2521 0 1.02962 -0.05154 0
2522 0 0.0595194 1.18892 0
2523 */
2524
2525 THD_3tta_to_3mni(&x, &y, &z); /* xform tt to mni space - results in lpi order */
2526 *xout = -x; *yout = -y; *zout = z; /* put coords back in RAI */
2527
2528 return(0);
2529 }
2530
2531
2532
2533 /* apply Brett transform to transform from MNI (MNI to Talairach)*/
2534 /* special case of 2-piece transform and maybe the only one we will ever use */
2535 int
apply_xform_brett_mni2tt(float x,float y,float z,float * xout,float * yout,float * zout)2536 apply_xform_brett_mni2tt(float x, float y, float z, \
2537 float *xout, float *yout, float *zout)
2538 {
2539 x = -x; y = -y; /* put coords in lpi from RAI first */
2540 THD_3mni_to_3tta(&x, &y, &z); /* xform mni to tt space - results in RAI */
2541
2542 *xout = x; *yout = y; *zout = z;
2543 return(0);
2544 }
2545
2546
2547 /* read xform from NIML attributes into xform_struct */
atlas_read_xform(NI_element * nel,ATLAS_XFORM * atlas_xf)2548 int atlas_read_xform(NI_element *nel, ATLAS_XFORM *atlas_xf)
2549 {
2550 float dist;
2551 int i;
2552 char *sptr;
2553
2554 if(wami_verb() > 2) {
2555 INFO_message("xform_name %s", NI_get_attribute(nel, "xform_name"));
2556 INFO_message("xform_type %s", NI_get_attribute(nel, "xform_type"));
2557 INFO_message("xform source %s", NI_get_attribute(nel, "source"));
2558 INFO_message("xform dest %s", NI_get_attribute(nel, "dest"));
2559 INFO_message("xform number of elements %d", nel->vec_num);
2560 INFO_message("xform post %s", NI_get_attribute(nel, "post"));
2561 INFO_message("xform coord_order %s", NI_get_attribute(nel, "coord_order"));
2562 }
2563 atlas_xf->xform_type = nifti_strdup(NI_get_attribute(nel, "xform_type"));
2564 atlas_xf->xform_name = nifti_strdup(NI_get_attribute(nel, "xform_name"));
2565 atlas_xf->source = nifti_strdup(NI_get_attribute(nel, "source"));
2566 atlas_xf->dest = nifti_strdup(NI_get_attribute(nel, "dest"));
2567 if(NI_get_attribute(nel, "distance")){
2568 dist = atof( NI_get_attribute(nel, "distance"));
2569 if(dist<=0) {
2570 WARNING_message("Distance less than or equal to 0 reset to 1");
2571 dist = 1;
2572 }
2573 }
2574 else
2575 dist = 1;
2576 atlas_xf->dist = dist;
2577
2578 sptr = NI_get_attribute(nel, "post");
2579 if(sptr){
2580 atlas_xf->post = atoi(sptr);
2581 }
2582 else
2583 atlas_xf->post = 0; /*assume pre-xform (used for 2 and 12 part xforms */
2584
2585 sptr = NI_get_attribute(nel, "coord_order");
2586 if(sptr){
2587 atlas_xf->coord_order = nifti_strdup(sptr);
2588 }
2589 else
2590 atlas_xf->coord_order = nifti_strdup("rai");
2591
2592 if((atlas_xf->xform_type == NULL) || (atlas_xf->source == NULL) ||
2593 (atlas_xf->dest == NULL) || (atlas_xf->xform_name==NULL) ||
2594 (atlas_xf->coord_order == NULL)) {
2595 WARNING_message("Could not allocate transformation type string");
2596 return(1);
2597 }
2598
2599 atlas_xf->nelts = nel->vec_num;
2600 atlas_xf->inverse = 0;
2601
2602 atlas_xf->xform = calloc(nel->vec_num, sizeof(float));
2603 if(atlas_xf->xform == NULL) {
2604 WARNING_message("Could not allocate transformation");
2605 return(1);
2606 }
2607
2608 /* ZSS To Daniel: There are cases where nel->vec is NULL but vec_num == 1
2609 I don't think this should happen. That is what was
2610 causing the crash. Now I check for it.
2611 Under certain situations, I can generate plenty of
2612 those instances. For example:
2613 cd /Users/ziad/SUMA_Class_New/data/suma_demo/afni
2614 3dROIstats -nzmean -nomeanout \
2615 -mask lh.OccROIs.FULL.niml.dset v2s.lh.TS.FULL.niml.dset
2616 */
2617 if (wami_verb() && nel->vec_num && !nel->vec) {
2618 WARNING_message("Strange xform nel: Have vec_num=%d but NULL nel->vec",
2619 nel->vec_num);
2620 }
2621 for(i=0;i<nel->vec_num && nel->vec;i++){
2622 /* ZSS to Daniel: This memcpy is suspect.
2623 You are copying only the first value
2624 from each vector vec[i]. Assuming you
2625 allocate for vec_len*vec_num you can do:
2626 memcpy((char *)(atlas_xf->xform)+(i*sizeof(float)*vec_len),
2627 nel->vec[i], sizeof(float)*vec_len);
2628 */
2629 memcpy((char *)(atlas_xf->xform)+i*sizeof(float),
2630 nel->vec[i], sizeof(float));
2631 }
2632 if(wami_verb() > 2)
2633 print_xform(atlas_xf);
2634 return(0);
2635 }
2636
2637 /* read template info from NIML attributes into template structure */
atlas_read_template(NI_element * nel,ATLAS_TEMPLATE * atlas_tpl)2638 int atlas_read_template(NI_element *nel, ATLAS_TEMPLATE *atlas_tpl)
2639 {
2640 char *s;
2641
2642 if(wami_verb() > 2) {
2643 INFO_message("template_name %s", NI_get_attribute(nel, "template_name"));
2644 INFO_message("templ_space %s", NI_get_attribute(nel, "template_space"));
2645 }
2646 atlas_tpl->comment = NULL;
2647 atlas_tpl->description = NULL;
2648
2649 atlas_tpl->template =
2650 nifti_strdup(NI_get_attribute(nel, "template_name"));
2651 atlas_tpl->space =
2652 nifti_strdup(NI_get_attribute(nel, "template_space"));
2653
2654 if ((s=NI_get_attribute(nel,"comment")))
2655 atlas_tpl->comment = nifti_strdup(s);
2656 if ((s=NI_get_attribute(nel,"description")))
2657 atlas_tpl->description = nifti_strdup(s);
2658
2659 if((atlas_tpl->template == NULL) || (atlas_tpl->space == NULL)) {
2660 WARNING_message("Could not get template strings");
2661 return(1);
2662 }
2663
2664 return(0);
2665 }
2666
2667 /* read atlas info from NIML attributes into atlas structure */
atlas_read_atlas(NI_element * nel,ATLAS * atlas,char * parentdir)2668 int atlas_read_atlas(NI_element *nel, ATLAS *atlas, char *parentdir)
2669 {
2670 char *s;
2671
2672 if (ATL_DSET(atlas) || atlas->adh || ATL_NAME(atlas)) {
2673 ERROR_message("Unclean atlas structure.");
2674 return(1);
2675 }
2676 if(wami_verb() > 2) {
2677 INFO_message("atlas_name %s", NI_get_attribute(nel, "atlas_name"));
2678 INFO_message("atlas_space %s", NI_get_attribute(nel, "template_space"));
2679 }
2680
2681 /* initialize the atlas fields */
2682 atlas->name = NULL;
2683 atlas->dset_name = NULL;
2684 atlas->comment = NULL;
2685 atlas->description = NULL;
2686 atlas->orient = NULL; /* coordinate order for web queries (Elsevier) */
2687 atlas->atlas_type = NULL; /* atlas can be "web"/ NULL */
2688 atlas->supp_web_info = NULL; /* supplemental info available at website */
2689 atlas->supp_web_type = NULL; /* atlas info may be .pdf, .html, ...*/
2690 atlas->supp_conn_info = NULL; /* supplemental connection info available at website */
2691 atlas->supp_conn_type = NULL; /* atlas info may be .pdf, .html, ...*/
2692 atlas->atlas_found = 0; /* flag for dataset available, -1 not found,
2693 1 found, 0 init value */
2694
2695 if ((s=NI_get_attribute(nel, "dset_name"))) {
2696 atlas->dset_name = NULL;
2697 if (!THD_is_prefix_ondisk(s, 0) &&
2698 parentdir && !THD_filehaspath(s)) {
2699 char *ss=(char *)calloc(strlen(parentdir)+strlen(s)+2,sizeof(char*));
2700 sprintf(ss,"%s/%s",parentdir,s);
2701 if (THD_is_prefix_ondisk(ss, 0))
2702 atlas->dset_name = nifti_strdup(ss);
2703 free(ss); ss=NULL;
2704 }
2705 if (!atlas->dset_name) atlas->dset_name = nifti_strdup(s);
2706 }
2707 if ((s=NI_get_attribute(nel, "template_space")))
2708 atlas->space = nifti_strdup(s);
2709 if ((s=NI_get_attribute(nel,"atlas_name")))
2710 atlas->name = nifti_strdup(s);
2711 if ((s=NI_get_attribute(nel,"description")))
2712 atlas->description = nifti_strdup(s);
2713 if ((s=NI_get_attribute(nel,"comment")))
2714 atlas->comment = nifti_strdup(s);
2715 if ((s=NI_get_attribute(nel,"orient")))
2716 atlas->orient = nifti_strdup(s);
2717 if ((s=NI_get_attribute(nel,"type")))
2718 atlas->atlas_type = nifti_strdup(s);
2719 if ((s=NI_get_attribute(nel,"supp_web_info")))
2720 atlas->supp_web_info = nifti_strdup(s);
2721 if ((s=NI_get_attribute(nel,"supp_web_type")))
2722 atlas->supp_web_type = nifti_strdup(s);
2723 if ((s=NI_get_attribute(nel,"supp_conn_info")))
2724 atlas->supp_conn_info = nifti_strdup(s);
2725 if ((s=NI_get_attribute(nel,"supp_conn_type")))
2726 atlas->supp_conn_type = nifti_strdup(s);
2727
2728 if((atlas->dset_name == NULL) || (atlas->space == NULL)) {
2729 WARNING_message("bad atlas nel");
2730 return(1);
2731 }
2732
2733 atlas->adh = NULL;
2734
2735 return(0);
2736 }
2737
2738 /* duplicate atlas info from one atlas structure to another */
2739 /* pointers are copied, no new allocation */
atlas_dup_atlas(ATLAS * srcatlas,ATLAS * destatlas)2740 int atlas_dup_atlas(ATLAS *srcatlas, ATLAS *destatlas)
2741 {
2742 /* initialize the atlas fields */
2743 destatlas->dset_name = srcatlas->dset_name;
2744 destatlas->space = srcatlas->space;
2745 destatlas->name = srcatlas->name;
2746 destatlas->description = srcatlas->description;
2747 destatlas->comment = srcatlas->comment;
2748 destatlas->atlas_found = srcatlas->atlas_found;
2749 destatlas->orient = srcatlas->orient;
2750 destatlas->atlas_type = srcatlas->atlas_type;
2751 destatlas->supp_web_info = srcatlas->supp_web_info;
2752 destatlas->supp_web_type = srcatlas->supp_web_type;
2753 destatlas->supp_conn_info = srcatlas->supp_conn_info;
2754 destatlas->supp_conn_type = srcatlas->supp_conn_type;
2755
2756 destatlas->adh = srcatlas->adh;
2757
2758 return(0);
2759 }
2760
2761 /* read template space info from NIML attributes into template space structure */
atlas_read_atlas_space(NI_element * nel,ATLAS_SPACE * at_space)2762 int atlas_read_atlas_space(NI_element *nel, ATLAS_SPACE *at_space)
2763 {
2764 if(wami_verb() > 2) {
2765 INFO_message("space_name %s", NI_get_attribute(nel, "space_name"));
2766 INFO_message("generic_space %s", NI_get_attribute(nel, "generic_space"));
2767 }
2768
2769 at_space->atlas_space = nifti_strdup(NI_get_attribute(nel, "space_name"));
2770 at_space->generic_space = nifti_strdup(NI_get_attribute(nel, "generic_space"));
2771
2772 if((at_space->atlas_space == NULL) || (at_space->generic_space == NULL)) {
2773 WARNING_message("Could not allocate template space strings");
2774 return(1);
2775 }
2776
2777 return(0);
2778 }
2779
2780
2781
NI_find_next_element(NI_stream ns,char * name)2782 NI_element *NI_find_next_element(NI_stream ns, char *name)
2783 {
2784 NI_element *nel;
2785
2786 nel = (NI_element *) 1;
2787 while(nel) {
2788 nel = NI_read_element(ns, 100);
2789 if(nel) {
2790 if (wami_verb()>2) fprintf(stderr,"nel name %s\n", nel->name);
2791 if (nel->type == NI_ELEMENT_TYPE) {
2792 if(strcmp(name, nel->name) == 0) {
2793 if (wami_verb()>2) fprintf(stderr, "name matches \n");
2794 return(nel);
2795 }
2796 }
2797 }
2798 }
2799 return(NULL);
2800 }
2801
2802 /* convert some of the atlas lists that are hard-coded in AFNI to NIML tables */
AFNI_atlas_list_to_niml()2803 void AFNI_atlas_list_to_niml()
2804 {
2805 ATLAS_POINT_LIST *temp_apl;
2806
2807 INFO_message("This is a debugging function. Not for regular use.");
2808
2809 /* if(wami_verb() > 1)
2810 INFO_message("Converting TTO_list_HARD to atlas_point_list");
2811 temp_apl = AFNI_atlas_list_to_atlas_point_list(TTO_list_HARD, TTO_COUNT_HARD);
2812 adjust_atlas_point_list(temp_apl, "Left", 200);
2813
2814 if(wami_verb() > 1){
2815 print_atlas_point_list(temp_apl);
2816 INFO_message("NIMLizing TTO_list_HARD");
2817 }
2818 atlas_list_to_niml(temp_apl, "TT_atlas.niml");
2819 if(wami_verb() > 1)
2820 INFO_message("Freeing the atlas_point_list");
2821 free_atlas_point_list(temp_apl);
2822 */
2823 if(wami_verb() > 2)
2824 INFO_message("Converting CA_EZ_list_HARD to atlas_point_list");
2825 temp_apl = AFNI_atlas_list_to_atlas_point_list(CA_EZ_list_HARD, CA_EZ_COUNT_HARD);
2826
2827 if(wami_verb() > 1){
2828 print_atlas_point_list(temp_apl);
2829 INFO_message("NIMLizing CA_EZ_list_HARD");
2830 }
2831 atlas_list_to_niml(temp_apl, "CA_EZ_atlas.niml");
2832 if(wami_verb() > 1)
2833 INFO_message("Freeing the atlas_point_list");
2834 free_atlas_point_list(temp_apl);
2835 /*
2836 if(wami_verb() > 1)
2837 INFO_message("Converting ML_EZ_list_HARD to atlas_point_list");
2838 temp_apl = AFNI_atlas_list_to_atlas_point_list(ML_EZ_list_HARD, ML_EZ_COUNT_HARD);
2839
2840 if(wami_verb() > 1){
2841 print_atlas_point_list(temp_apl);
2842 INFO_message("NIMLizing ML_EZ_list_HARD");
2843 }
2844 atlas_list_to_niml(temp_apl, "ML_EZ_atlas.niml");
2845 if(wami_verb() > 1)
2846 INFO_message("Freeing the atlas_point_list");
2847 free_atlas_point_list(temp_apl);
2848
2849 if(wami_verb() > 1)
2850 INFO_message("Converting LR_EZ_list_HARD to atlas_point_list");
2851 temp_apl = AFNI_atlas_list_to_atlas_point_list(LR_EZ_list_HARD, LR_EZ_COUNT_HARD);
2852
2853 if(wami_verb() > 1){
2854 print_atlas_point_list(temp_apl);
2855 INFO_message("NIMLizing LR_EZ_list_HARD");
2856 }
2857 atlas_list_to_niml(temp_apl, "LR_EZ_atlas.niml");
2858 if(wami_verb() > 1)
2859 INFO_message("Freeing the atlas_point_list");
2860 free_atlas_point_list(temp_apl);
2861
2862 */
2863 }
2864
2865 #define TRIM_STRING(lbuf, ch){\
2866 int kk; \
2867 for( kk=strlen(lbuf)-1 ; kk > 0 && lbuf[kk] == ch ; kk-- ) \
2868 lbuf[kk] = '\0' ; /* trim trailing .'s */ \
2869 }
2870 /* this function was used originally to convert an old atlas format
2871 * that was hard coded here to the newer one used in AFNI header
2872 * attributes. This is likely not to be used again.
2873 */
AFNI_atlas_list_to_atlas_point_list(ATLAS_POINT * afni_at_pts,int npts)2874 static ATLAS_POINT_LIST * AFNI_atlas_list_to_atlas_point_list(
2875 ATLAS_POINT *afni_at_pts, int npts)
2876 {
2877 ATLAS_POINT *temp_atp;
2878 ATLAS_POINT_LIST *apl;
2879 int i;
2880
2881 ENTRY("AFNI_atlas_list_to_atlas_point_list");
2882 apl = (ATLAS_POINT_LIST *) calloc(1,sizeof(ATLAS_POINT_LIST));
2883 apl->n_points = npts;
2884 apl->at_point = (ATLAS_POINT *) calloc(npts,sizeof(ATLAS_POINT));
2885 for(i=0;i<npts;i++){
2886 /* copy AFNI point value to new atlas point list structure */
2887 temp_atp = apl->at_point+i;
2888 temp_atp->tdval = afni_at_pts[i].tdval;
2889 temp_atp->tdlev = afni_at_pts[i].tdlev;
2890 temp_atp->okey = afni_at_pts[i].okey;
2891 temp_atp->xx = afni_at_pts[i].xx;
2892 temp_atp->yy = afni_at_pts[i].yy;
2893 temp_atp->zz = afni_at_pts[i].zz;
2894 NI_strncpy(temp_atp->name,afni_at_pts[i].name,ATLAS_CMAX);
2895 TRIM_STRING(temp_atp->name, '.');
2896 NI_strncpy(temp_atp->sblabel,afni_at_pts[i].sblabel,ATLAS_CMAX);
2897 TRIM_STRING(temp_atp->sblabel, '.');
2898 NI_strncpy(temp_atp->longname, "", ATLAS_CMAX);
2899 if(wami_verb() > 1){
2900 INFO_message("atlas_point %d %s\n", afni_at_pts[i].tdval,
2901 afni_at_pts[i].name);
2902 INFO_message("atlas_point %d %s temp\n", temp_atp->tdval,
2903 temp_atp->name);
2904 }
2905 }
2906 RETURN(apl);
2907 }
2908
label_table_to_atlas_point_list(Dtable * dtbl)2909 ATLAS_POINT_LIST * label_table_to_atlas_point_list(Dtable *dtbl)
2910 {
2911 ATLAS_POINT *temp_atp;
2912 ATLAS_POINT_LIST *apl;
2913 int i, nn;
2914 char **la , **lb;
2915
2916 ENTRY("label_table_to_atlas_point_list");
2917
2918 nn = listize_Dtable( dtbl , &la , &lb ) ;
2919 if( nn == 0 || la == NULL || lb == NULL ) RETURN (NULL) ;
2920
2921 apl = (ATLAS_POINT_LIST *) calloc(1,sizeof(ATLAS_POINT_LIST));
2922 apl->n_points = nn;
2923 apl->at_point = (ATLAS_POINT *) calloc(nn,sizeof(ATLAS_POINT));
2924 for(i=0;i<nn;i++){
2925 temp_atp = apl->at_point+i;
2926 temp_atp->tdval = (int)strtol(la[i],NULL,10);
2927 temp_atp->tdlev = 0;
2928 temp_atp->okey = (int)strtol(la[i],NULL,10);
2929 temp_atp->xx = 0.0;
2930 temp_atp->yy = 0.0;
2931 temp_atp->zz = 0.0;
2932 NI_strncpy(temp_atp->name,lb[i],ATLAS_CMAX);
2933 TRIM_STRING(temp_atp->name, '.');
2934 NI_strncpy(temp_atp->sblabel,lb[i],ATLAS_CMAX);
2935 TRIM_STRING(temp_atp->sblabel, '.');
2936 NI_strncpy(temp_atp->longname,"",ATLAS_CMAX);
2937 if(wami_verb() > 1){
2938 INFO_message("Dtable %d %s\n", (int)strtol(la[i],NULL,10),
2939 lb[i]);
2940 INFO_message("atlas_point %d %s temp\n", temp_atp->tdval,
2941 temp_atp->name);
2942 }
2943 }
2944 RETURN(apl);
2945 }
2946
2947 /* convert atlas list to a NIML table in a text file */
2948 void
atlas_list_to_niml(ATLAS_POINT_LIST * atp,char * atlas_file)2949 atlas_list_to_niml(ATLAS_POINT_LIST *atp, char *atlas_file)
2950 {
2951 int i;
2952 /* char *atlas_pt_niml;*/
2953 char filestr[ATLAS_CMAX];
2954
2955 ATLAS_POINT *at_pt;
2956 NI_stream atlas_niml;
2957 NI_group *ngr;
2958 NI_element *nel=NULL;
2959
2960 ENTRY("atlas_list_to_niml");
2961
2962 if(wami_verb() > 1)
2963 INFO_message("opening %s", atlas_file);
2964 if (atlas_file == NULL) {
2965 sprintf(filestr, "stdout:");
2966 } else {
2967 sprintf(filestr, "file:%s", atlas_file);
2968 }
2969 atlas_niml = NI_stream_open(filestr,"w");
2970 if(atlas_niml==NULL){
2971 WARNING_message("Could not open atlas file for NIML output %s");
2972 EXRETURN;
2973 }
2974 ngr = NI_new_group_element();
2975 NI_rename_group(ngr, "atlas_point_list");
2976 /* get each segmented region - the "atlas point" into
2977 a NIML formatted string */
2978 for(i=0;i<atp->n_points;i++) {
2979 at_pt = atp->at_point+i;
2980 nel = atlas_point_to_niml_element(at_pt);
2981 NI_add_to_group( ngr, nel);
2982 /* atlas_pt_niml = atlas_point_to_niml_string(at_pt);*/
2983 /* write the NIML string to the NIML file */
2984 /* if(write_niml_string(atlas_pt_niml, atlas_niml)) {
2985 WARNING_message("Could not write atlas point to NIML file");
2986 NI_stream_close(atlas_niml);
2987 EXRETURN;
2988 }
2989 */
2990 }
2991
2992 /* write the whole group out at once to a NIML file */
2993 if(NI_write_element( atlas_niml, ngr , NI_TEXT_MODE&NI_HEADERSHARP_FLAG)<0) {
2994 WARNING_message("Could not write atlas point list to NIML file");
2995 }
2996
2997 NI_free_element(ngr) ;
2998 NI_stream_close(atlas_niml);
2999 EXRETURN;
3000 }
3001
3002 /* adjust the point values in an atlas list by a fixed amount if the name */
3003 /* of the structure starts with the specified match string - case insensitive*/
adjust_atlas_point_list(ATLAS_POINT_LIST * atp,char * match_str,float addval)3004 static void adjust_atlas_point_list(ATLAS_POINT_LIST *atp, char *match_str, float addval)
3005 {
3006 /* this function was written for adjusting the TT_daemon values to be unique
3007 previously left and right regions shared the same value */
3008 int i;
3009 ATLAS_POINT *at_pt;
3010
3011 ENTRY("adjust_atlas_point_list");
3012
3013 for(i=0;i<atp->n_points;i++) {
3014 at_pt = atp->at_point+i;
3015 if((strstr(at_pt->name, match_str))){
3016 at_pt->tdval = at_pt->tdval + addval;
3017 }
3018 }
3019
3020 EXRETURN;
3021 }
3022
3023
3024
3025 #if 0
3026 /* write a niml encoded string to a specified stream */
3027 static int
3028 write_niml_string(char *nimlstring, NI_stream ns)
3029 {
3030 NI_element *nel=NULL;
3031 ENTRY("write_niml_string");
3032 nel = (NI_element *)NI_read_element_fromstring(nimlstring);
3033 /* write the niml element into a string - note this is written in text mode
3034 (ascii encoded) */
3035 if(NI_write_element( ns, nel , NI_TEXT_MODE&NI_HEADERSHARP_FLAG)<0)
3036 RETURN(-1);
3037 RETURN(0);
3038 }
3039 #endif
3040
3041 #if 0
3042 just for reference here
3043 typedef struct {
3044 /* tdval and tdlev stand for "Talairach Daemon" value and level */
3045 /* these are kept for historical purposes */
3046 /* perhaps one day making an unusally boring PBS special */
3047 short tdval; /* Leave this one to be the very first element */
3048 char name[ATLAS_CMAX] ; /* Leave this one to be the second element */
3049 short xx,yy,zz,tdlev,okey ; /* xx,yy,zz - RAI position of region */
3050 /* tdlev = unknown, gyrus or area code */
3051 /* okey = original value in atlas */
3052 /* this value was converted for TT daemon */
3053 /* atlas values because left and right */
3054 /* ROIs shared the same value */
3055 char sblabel[ATLAS_CMAX];
3056 } ATLAS_POINT ;
3057 #endif
3058
3059 /* return niml element from atlas point structure */
3060 NI_element *
atlas_point_to_niml_element(ATLAS_POINT * at_pt)3061 atlas_point_to_niml_element(ATLAS_POINT *at_pt)
3062 {
3063 float cog[3];
3064 NI_element *nel=NULL;
3065 short okey, tdval;
3066
3067 ENTRY("atlas_point_to_niml_element");
3068
3069 /* create niml element that corresponds to atlas point structure */
3070 nel = NI_new_data_element("ATLAS_POINT",0);
3071 NI_set_attribute( nel,"data_type", "atlas_point");
3072 NI_set_attribute(nel, "STRUCT", at_pt->name);
3073
3074 /* original key value because I will be modifying the original values -
3075 need to put this into TTO_list */
3076 /* TT atlas includes the same values for left and right structures,
3077 but I'll be adding 200 to these */
3078 tdval = at_pt->tdval;
3079 NI_SET_INT(nel, "VAL", tdval);
3080 okey = at_pt->okey;
3081 if(okey == -999) okey = tdval; /* get original "key" value from tdval */
3082 NI_SET_INT(nel, "OKEY", okey);
3083
3084 /* gyrus or area code */
3085 NI_SET_INT(nel, "GYoAR", at_pt->tdlev);
3086
3087 /* center of mass/gravity coordinates */
3088 cog[0] = at_pt->xx; cog[1] = at_pt->yy; cog[2] = at_pt->zz;
3089 NI_SET_FLOATv(nel, "COG", cog, 3);
3090
3091 /* associated dataset file that gives the original name of the Analyze file
3092 suggesting the name of the sub-brick in the atlas where
3093 it is finally stored */
3094 if (strcmp(at_pt->sblabel,"")!=0)
3095 NI_set_attribute(nel, "SB_LABEL", at_pt->sblabel);
3096 if (strcmp(at_pt->longname,"")!=0)
3097 NI_set_attribute(nel, "LONGNAME", at_pt->longname);
3098
3099 RETURN(nel);
3100 }
3101
3102 /* return ascii encoded niml string from atlas point structure */
3103 char *
atlas_point_to_niml_string(ATLAS_POINT * at_pt)3104 atlas_point_to_niml_string(ATLAS_POINT *at_pt)
3105 {
3106 NI_element *nel=NULL;
3107 NI_stream ns;
3108 char *encstr;
3109 nel = atlas_point_to_niml_element(at_pt);
3110 /* create a niml stream to write into */
3111 ns = NI_stream_open("str:", "w");
3112
3113 if(ns==NULL)
3114 RETURN(NULL);
3115 /* write the niml element into a string -
3116 note this is written in text mode (ascii encoded) */
3117 NI_write_element( ns , nel , NI_TEXT_MODE&NI_HEADERSHARP_FLAG) ;
3118
3119 /* copy the string to a more permanent string and
3120 then can close stream and string */
3121 encstr = SUMA_copy_string( NI_stream_getbuf(ns));
3122
3123 /* close stream - freeing elements */
3124 NI_stream_close( ns );
3125
3126 RETURN(encstr);
3127 }
3128
3129 /* reading functions */
3130
3131 /* convert a NIML table from a file to an atlas list structure */
3132 static void
niml_to_atlas_list(ATLAS_POINT_LIST * atp,char * atlas_file)3133 niml_to_atlas_list(ATLAS_POINT_LIST *atp, char *atlas_file)
3134 {
3135 int more_niml, count, tdlev;
3136 char filestr[ATLAS_CMAX];
3137 char *temp_str;
3138 ATLAS_POINT *at_pt;
3139 NI_stream atlas_niml;
3140 NI_element *nel=NULL;
3141
3142 float cog[3];
3143 char *encstr=NULL;
3144 short okey;
3145
3146 ENTRY("niml_to_atlas_list");
3147
3148 if(wami_verb() > 1)
3149 INFO_message("opening atlas niml file %s",atlas_file);
3150 sprintf(filestr, "file:%s", atlas_file);
3151 atlas_niml = NI_stream_open(filestr,"r");
3152
3153 if(atlas_niml==NULL){
3154 WARNING_message("Could not open atlas file for NIML output %s");
3155 EXRETURN;
3156 }
3157
3158 more_niml = 1;
3159 count = 0;
3160 /* get each segmented region - the "atlas point" from
3161 a NIML formatted string */
3162 while(more_niml) {
3163 /* read the NIML string from the NIML file */
3164 nel = (NI_element *)NI_read_element_fromstring(encstr);
3165 if(!nel) {
3166 NI_stream_close(atlas_niml);
3167 EXRETURN;
3168 }
3169 NI_GET_INT(nel, "OKEY", okey);
3170 NI_GET_INT(nel, "GYoAR", tdlev);
3171 NI_GET_FLOATv(nel, "COG", cog, 3, 0);
3172
3173 temp_str = NI_get_attribute(nel, "STRUCT");
3174 at_pt = &atp->at_point[count];
3175 NI_strncpy(at_pt->name,temp_str,ATLAS_CMAX);
3176
3177 temp_str = NI_get_attribute(nel, "SB_LABEL");
3178 if(temp_str==NULL)
3179 NI_strncpy(at_pt->sblabel,"",ATLAS_CMAX);
3180 else
3181 NI_strncpy(at_pt->sblabel,temp_str,ATLAS_CMAX);
3182
3183 temp_str = NI_get_attribute(nel, "LONGNAME");
3184 if(temp_str==NULL)
3185 NI_strncpy(at_pt->longname,"",ATLAS_CMAX);
3186 else
3187 NI_strncpy(at_pt->longname,temp_str,ATLAS_CMAX);
3188
3189 at_pt->tdval = okey;
3190 at_pt->xx = cog[0];
3191 at_pt->yy = cog[1];
3192 at_pt->zz = cog[2];
3193 if(wami_verb() > 1)
3194 fprintf(stderr,"Decoded:\n"
3195 "STRUCT: %s\n"
3196 "OKEY: %d\n"
3197 "COG: %.3f %.3f %.3f\n"
3198 "SB_LABEL: %s\n"
3199 , NI_get_attribute(nel, "STRUCT"),
3200 okey, cog[0], cog[1], cog[2],
3201 NI_get_attribute(nel, "SB_LABEL"));
3202 count++;
3203 }
3204 NI_stream_close(atlas_niml);
3205 EXRETURN;
3206 }
3207
3208 #if 0
3209 static int
3210 read_niml_string() /* Now pretend we got encstr from our label table, by using the key
3211 {
3212 okey, or okey+256 */
3213 EXRETURN;
3214 }
3215 #endif
3216
3217 /* get default jump to space, usually "jump to MNI" for AFNI gui */
get_jump_space()3218 char *get_jump_space()
3219 {
3220 char *spacename;
3221
3222 if(jumpspace_name)
3223 return(jumpspace_name);
3224
3225 spacename = getenv("AFNI_JUMPTO_SPACE") ;
3226 if(spacename && (strlen(spacename)>0) && (strlen(spacename)<110)){
3227 jumpspace_name = nifti_strdup(spacename);
3228 }
3229 else
3230 jumpspace_name = nifti_strdup("MNI");
3231
3232 return(jumpspace_name);
3233 }
3234
3235 /* get default jump to space, usually "jump to MNI" */
set_jump_space(char * spacename)3236 void set_jump_space(char *spacename)
3237 {
3238 if(spacename && (strlen(spacename)>0) && (strlen(spacename)<110)){
3239 if(jumpspace_name) free(jumpspace_name);
3240 jumpspace_name = nifti_strdup(spacename);
3241 }
3242 }
3243
3244
3245