1 /*! Dealings with volume data */
2 #include "SUMA_suma.h"
3
4 extern SUMA_CommonFields *SUMAg_CF;
5
SUMA_WarpTypeName(SUMA_WARP_TYPES wt)6 const char *SUMA_WarpTypeName(SUMA_WARP_TYPES wt)
7 {
8 switch (wt) {
9 case BAD_WARP:
10 return ("Bad Warp");
11 case NO_WARP:
12 return ("No Warp");
13 case ROTATE_WARP:
14 return ("3drotate Warp");
15 case VOLREG_WARP:
16 return ("3dvolreg Warp");
17 case ALLINEATE_WARP:
18 return ("3dAllineate Warp");
19 case WARPDRIVE_WARP:
20 return ("3dWarpdrive Warp");
21 case TAGALIGN_WARP:
22 return ("3dTagalign Warp");
23 case N_WARP_TYPES:
24 return ("Number of Warp Types");
25 default:
26 return ("You're a warped type");
27 }
28 }
29
30 /*!
31 \brief Find the closest node on the surface to a voxel in the volume (as defined in vp)
32 \param SO (SUMA_SurfaceObject *)
33 SO->NodeList is expected to be in dicomm (RAI) units
34 \param vp (SUMA_VOLPAR *) Volume grid of Nvox voxels.
35 \param closest_node (int *) Vector of Nvox elements. To contain the nodes closest to each voxel
36 \param closest_dist (float *) Vector of Nvox elements.
37 To contain the distance between node and voxel centroid
38 Pass NULL if you do not care for this info.
39 \param vox_mask (byte *) Mask of voxels to process (NULL if you want to process all)
40 \param verb (int): 0 be quiet
41 > 1 talk
42
43 */
SUMA_ClosestNodeToVoxels(SUMA_SurfaceObject * SO,SUMA_VOLPAR * vp,int * closest_node,float * closest_dist,byte * vox_mask,int verb)44 int SUMA_ClosestNodeToVoxels(SUMA_SurfaceObject *SO, SUMA_VOLPAR *vp,
45 int *closest_node, float *closest_dist,
46 byte *vox_mask, int verb)
47 {
48 static char FuncName[]={"SUMA_ClosestNodeToVoxels"};
49 float *p=NULL;
50 float d, dxyz;
51 int i, j, k, n, nij, ijk, cnt = 0;
52 THD_fvec3 fv , iv;
53
54 SUMA_ENTRY;
55
56 if (!SO || !vp || !closest_node) {
57 SUMA_S_Err("NULL input");
58 SUMA_RETURN(0);
59 }
60
61 if (verb) {
62 fprintf(SUMA_STDERR, "%s: Have %d nodes in surface\n"
63 "%dx%dx%d (%d) voxels in volume.\n",
64 FuncName, SO->N_Node, vp->nx, vp->ny, vp->nz,
65 vp->nx * vp->ny * vp->nz);
66 }
67 /* Now for each voxel, find the closest node (SLOW implementation) */
68 cnt = 0;
69 nij = vp->nx*vp->ny;
70 for (i=0;i<vp->nx; ++i) {
71 for (j=0;j<vp->ny; ++j) {
72 for (k=0;k<vp->nz; ++k) {
73 ijk = SUMA_3D_2_1D_index(i,j,k,vp->nx,nij);
74 /* fprintf(SUMA_STDERR," %4d %4d %4d, ijk: %d\n", i, j, k, ijk); */
75 closest_node[ijk] = -1;
76 if (closest_dist) closest_dist[ijk] = -1.0;
77 if (!vox_mask || (vox_mask && vox_mask[ijk])) {
78 iv.xyz[0] = (float)i; iv.xyz[1] = (float)j; iv.xyz[2] = (float)k;
79 fv = SUMA_THD_3dfind_to_3dmm_vp(vp, iv);
80 iv = SUMA_THD_3dmm_to_dicomm(vp->xxorient, vp->yyorient,
81 vp->zzorient, fv);
82 #if 0 /* macro not tested here so keeping older code below */
83 SUMA_CLOSEST_NODE(SO, iv.xyz, closest_node[ijk], d);
84 if (closest_dist) closest_dist[ijk] = (float)d;
85 #else
86 dxyz = 1023734552736672366372.0;
87 for (n=0; n<SO->N_Node; ++n) {
88 p = &(SO->NodeList[SO->NodeDim*n]);
89 SUMA_SEG_LENGTH_SQ(p, iv.xyz, d);
90 if (d < dxyz) {
91 dxyz = d; closest_node[ijk] = n;
92 if (closest_dist) closest_dist[ijk] = (float)d;
93 }
94 }
95 #endif
96 if (closest_dist) {
97 if (closest_dist[ijk] >= 0.0f)
98 closest_dist[ijk] = (float)sqrt(closest_dist[ijk]); }
99 if (verb) {
100 ++cnt;
101 if (!(cnt % 1000)) {
102 fprintf(SUMA_STDERR,". @ %4d %4d %4d (%3.2f%%)\n",
103 i, j, k,
104 (float)cnt/(float)(vp->nx * vp->ny * vp->nz)*100.0);
105 fflush(SUMA_STDERR);
106 }
107 }
108 }
109 }
110 }
111 }
112
113 SUMA_RETURN(1);
114 }
115
116
SUMA_THD_handedness(THD_3dim_dataset * dset)117 int SUMA_THD_handedness( THD_3dim_dataset * dset )
118 {
119 static char FuncName[]={"SUMA_THD_handedness"};
120
121 SUMA_ENTRY;
122
123 SUMA_RETURN(THD_handedness(dset));
124 }
125
126 /*!
127 see help for SUMA_AfniPrefix below
128 */
SUMA_AfniExistsView(int exists,char * view)129 SUMA_Boolean SUMA_AfniExistsView(int exists, char *view)
130 {
131 static char FuncName[]={"SUMA_AfniExistsView"};
132 SUMA_Boolean ans = NOPE;
133
134 SUMA_ENTRY;
135
136 if (!exists) SUMA_RETURN(ans);
137
138 if (strcmp(view,"+orig") == 0) {
139 if (exists == 1 || exists == 3 || exists == 5 || exists == 7) ans = YUP;
140 } else if (strcmp(view,"+acpc") == 0) {
141 if (exists == 2 || exists == 3 || exists == 6 || exists == 7) ans = YUP;
142 } else if (strcmp(view,"+tlrc") == 0) {
143 if (exists == 4 || exists == 5 || exists == 6 || exists == 7) ans = YUP;
144 }
145
146 SUMA_RETURN(ans);
147
148 }
149
SUMA_AfniView(char * nameorig,char * cview)150 SUMA_Boolean SUMA_AfniView (char *nameorig, char *cview)
151 {
152 static char FuncName[]={"SUMA_AfniView"};
153 char *tmp1 = NULL, *tmp2 = NULL;
154 SUMA_Boolean LocalHead = NOPE;
155
156 SUMA_ENTRY;
157
158 if (!nameorig) SUMA_RETURN(NOPE);
159 if (!cview) SUMA_RETURN(NOPE);
160
161 tmp1 = SUMA_Extension(nameorig, ".HEAD", YUP);
162 tmp2 = SUMA_Extension(tmp1, ".BRIK", YUP); SUMA_free(tmp1); tmp1 = NULL;
163 /* is there a dot ?*/
164 if (tmp2[strlen(tmp2)-1] == '.') tmp2[strlen(tmp2)-1] = '\0';
165
166 if (LocalHead) fprintf(SUMA_STDERR,"%s: Searching for view of %s\n", FuncName, tmp2);
167
168 /* view */
169 if (SUMA_isExtension(tmp2, "+orig")) {
170 sprintf(cview, "+orig");
171 } else if (SUMA_isExtension(tmp2, "+acpc")) {
172 sprintf(cview, "+acpc");
173 } else if (SUMA_isExtension(tmp2, "+tlrc")) {
174 sprintf(cview, "+tlrc");
175 } else {
176 cview[0]='\0';
177 }
178 SUMA_free(tmp2); tmp2 = NULL;
179
180 SUMA_RETURN(YUP);
181 }
182
SUMA_AfniExists(char * prefix,char * c2view)183 SUMA_Boolean SUMA_AfniExists(char *prefix, char *c2view)
184 {
185 static char FuncName[]={"SUMA_AfniExists"};
186 char *head=NULL;
187 SUMA_PARSED_NAME *Fname=NULL;
188 SUMA_Boolean ans = NOPE;
189 SUMA_Boolean LocalHead = NOPE;
190
191 SUMA_ENTRY;
192
193 ans = NOPE;
194
195 /* Jan 2015: Rely on ParseFname for all the work...*/
196 if (!(Fname = SUMA_ParseFname_eng (prefix, NULL,1))) {
197 SUMA_RETURN(ans);
198 }
199 head = SUMA_ModifyName(Fname->HeadName, "view",c2view, NULL);
200
201 SUMA_LH("Checking existence of %s\n", head);
202 if (SUMA_filexists(head)) { ans = YUP; }
203 SUMA_ifree(head); head = NULL;
204 SUMA_Free_Parsed_Name(Fname);
205
206 SUMA_RETURN(ans);
207 }
208 /*!
209 \brief Return AFNI's prefix, from a file name also checks for its validity
210 \param name (char *) dset name (can contain path)
211 \param view (char[4]) array to return view in it
212 \param path (char *) if any, make sure it is not duplicated in name!
213 \param exists (int *) function checks
214 for all three possible views
215 0: No afni dset with name/prefix
216 1: +orig
217 2: +acpc
218 3: +orig and +acpc
219 4: +tlrc
220 5: +orig and +tlrc
221 6: +acpc and +tlrc
222 7: +orig and +acpc and +tlrc
223 \sa SUMA_AfniExistsView
224
225 \return prefix (char *) dset prefix, free it with SUMA_free
226 \sa SUMA_AfniView
227 */
SUMA_AfniPrefix(char * nameorig,char * view,char * path,int * exists)228 char *SUMA_AfniPrefix(char *nameorig, char *view, char *path, int *exists)
229 {
230 static char FuncName[]={"SUMA_AfniPrefix"};
231 char *prfx = NULL, *name=NULL;
232 char cview[10];
233 int iview;
234 SUMA_PARSED_NAME *Fname=NULL;
235 SUMA_Boolean LocalHead = NOPE;
236
237 SUMA_ENTRY;
238
239 if (!nameorig) SUMA_RETURN(prfx);
240 if (exists) *exists = 0;
241
242 SUMA_LH("Working with %s\n", nameorig);
243
244 if (!path) name = SUMA_copy_string(nameorig);
245 else {
246 if (path[strlen(path)-1] == '/') {
247 name = SUMA_append_replace_string(path, nameorig, "", 0);
248 } else {
249 name = SUMA_append_replace_string(path, nameorig, "/", 0);
250 }
251 }
252
253 if (LocalHead) fprintf(SUMA_STDERR,"%s: name now %s\n", FuncName, name);
254
255 /* Jan 2015: Rely on ParseFname for all the work...*/
256 Fname = SUMA_ParseFname_eng (name, NULL,1);
257
258 /* view */
259 iview = -1;
260 if (!strcmp(Fname->View,"+orig")) {
261 iview = 0; sprintf(cview, "+orig");
262 } else if (!strcmp(Fname->View, "+acpc")) {
263 iview = 1; sprintf(cview, "+acpc");
264 } else if (!strcmp(Fname->View, "+tlrc")) {
265 iview = 2; sprintf(cview, "+tlrc");
266 } else {
267 cview[0]='\0';
268 }
269 prfx = SUMA_append_replace_string(Fname->Path,Fname->Prefix,"",0);
270
271 if (LocalHead) fprintf(SUMA_STDERR,"%s: Prefix %s\n", FuncName, prfx);
272
273 /* can't tell what view is, so can't test properly quite yet */
274 {
275 /* is name ok ? all I have to do is test with +orig */
276
277 if( !THD_filename_ok(Fname->HeadName) ){
278 SUMA_SL_Err("%s not a proper dset name!\n", Fname->HeadName) ;
279 SUMA_ifree(name); name = NULL;
280 SUMA_ifree(prfx); prfx = NULL;
281 SUMA_Free_Parsed_Name(Fname); Fname = NULL;
282 SUMA_RETURN(prfx);
283 }
284 }
285 if (exists) *exists = Fname->OnDisk;
286
287 if (exists) {
288 { /* look for any */
289 char *head=NULL, c2view[10];
290 SUMA_PARSED_NAME *Fname2=NULL;
291 iview = 0; *exists = 0;
292 SUMA_LH("Modifying %s\n", Fname->HeadName);
293 while (iview < 3) {
294 if (iview == 0) head =
295 SUMA_ModifyName(Fname->HeadName, "view","+orig", NULL);
296 else if (iview == 1) head =
297 SUMA_ModifyName(Fname->HeadName, "view","+acpc", NULL);
298 else if (iview == 2) head =
299 SUMA_ModifyName(Fname->HeadName, "view","+tlrc", NULL);
300 Fname2 = SUMA_ParseFname_eng (head, NULL,1);
301 SUMA_LH("Checking existence of %s\n", head);
302 if (Fname2->OnDisk) { *exists += (int)pow(2,iview); }
303 SUMA_free(head); head = NULL;
304 SUMA_Free_Parsed_Name(Fname2); Fname2 = NULL;
305 ++iview;
306 }
307 }
308 SUMA_LH("exist number is %d\n", *exists);
309 }
310
311
312 if (cview[0] != '\0') {
313 if (LocalHead) {
314 if (Fname->OnDisk) {
315 fprintf(SUMA_STDERR,"%s: dset with view %s does exist.\n",
316 FuncName, cview);
317 } else {
318 fprintf(SUMA_STDERR,"%s: dset with view %s does not exist.\n",
319 FuncName, cview);
320 }
321 }
322 }
323
324 if (view) {
325 sprintf(view,"%s", cview);
326 }
327 if (name) SUMA_free(name); name = NULL;
328 SUMA_Free_Parsed_Name(Fname); Fname = NULL;
329
330 SUMA_RETURN(prfx);
331 }
332
333 /*!
334 \brief A function to find the skin of a volume
335 \param dset (THD_3dim_dataset *) an AFNI volume
336 \param fvec (float *) (nx * ny * nz) data vector
337 \param thresh (double) consider only values in fvec > thresh
338 \param N_skin (int *) number of voxels that are skin
339 \return skin (byte *) (nx * ny * nz) vector containing 1 for skin voxels, 0 elsewhere.
340 */
SUMA_isSkin(THD_3dim_dataset * dset,float * fvec,double thresh,int * N_skin)341 byte * SUMA_isSkin(THD_3dim_dataset *dset, float *fvec, double thresh, int *N_skin)
342 {
343 static char FuncName[]={"SUMA_isSkin"};
344 byte *isskin=NULL;
345 int x, y, z, nx, ny, nz, i1D, nxy;
346
347 SUMA_ENTRY;
348
349 if (!dset || !fvec) {
350 SUMA_SL_Err("NULL input dset or fvec");
351 SUMA_RETURN(isskin);
352 }
353
354 nx = DSET_NX(dset);
355 ny = DSET_NY(dset);
356 nz = DSET_NZ(dset);
357 nxy = nx * ny;
358
359 isskin = (byte *) SUMA_calloc(nxy * nz, sizeof(byte));
360 if (!isskin) {
361 SUMA_SL_Crit("Failed to allocate");
362 SUMA_RETURN(NULL);
363 }
364
365 *N_skin = 0;
366 for (z=0; z<nz; ++z) {
367 for (y=0; y<ny; ++y) {
368 x = 0;
369 do { /* the upstroke */
370 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
371 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
372 ++x;
373 } while (x < nx && !isskin[i1D]);
374 x = nx - 1;
375 do { /* the downstroke */
376 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
377 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
378 --x;
379 } while (x >=0 && !isskin[i1D]);
380 } /* y */
381 } /* z */
382
383 for (z=0; z<nz; ++z) {
384 for (x=0; x<nx; ++x) {
385 y = 0;
386 do { /* the upstroke */
387 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
388 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
389 ++y;
390 } while (y < ny && !isskin[i1D]);
391 y = ny - 1;
392 do { /* the downstroke */
393 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
394 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
395 --y;
396 } while (y >=0 && !isskin[i1D]);
397 } /* x */
398 } /* z */
399
400 for (x=0; x<nx; ++x) {
401 for (y=0; y<ny; ++y) {
402 z = 0;
403 do { /* the upstroke */
404 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
405 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
406 ++z;
407 } while (z < nz && !isskin[i1D]);
408 z = nz - 1;
409 do { /* the downstroke */
410 i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
411 if (fvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
412 --z;
413 } while (z >=0 && !isskin[i1D]);
414 } /* y */
415 } /* x */
416
417 SUMA_RETURN(isskin);
418 }
419 /*!
420 Returns the attributes of the surface's volume parent
421 A volume parent is defined as:
422 The volume data set used to generate the surface. aka the master SPGR
423 \ret NULL for troubles
424 */
SUMA_Alloc_VolPar(void)425 SUMA_VOLPAR *SUMA_Alloc_VolPar (void)
426 {
427 static char FuncName[]={"SUMA_Alloc_VolPar"};
428 SUMA_VOLPAR *VP;
429
430 SUMA_ENTRY;
431
432 VP = (SUMA_VOLPAR *)SUMA_malloc(sizeof(SUMA_VOLPAR));
433 if (VP == NULL) {
434 fprintf(SUMA_STDERR,"Error SUMA_Alloc_VolPar: Failed to allocate for VolPar\n");
435 SUMA_RETURN (NULL);
436 }
437 VP->idcode_str = NULL;
438 VP->isanat = 1;
439 VP->nx = VP->ny = VP->nz = 0; /*!< number of voxels in the three dimensions */
440 VP->dx = VP->dy = VP->dz = 0.0; /*!< delta x, y, z in mm */
441 VP->xorg = VP->yorg = VP->zorg = 0.0; /*!< voxel origin in three dimensions */
442 VP->prefix = NULL; /*!< parent volume prefix */
443 VP->headname = NULL;
444 VP->filecode = NULL; /*!< parent volume prefix + view */
445 VP->dirname = NULL; /*!< parent volume directory name */
446 VP->vol_idcode_str = NULL; /*!< idcode string OF parent volume*/
447 VP->vol_idcode_date = NULL; /*!< idcode date */
448 VP->xxorient = VP->yyorient = VP->zzorient = 0; /*!< orientation of three dimensions*/
449 VP->CENTER_OLD = NULL; /*!< pointer to the named attribute (3x1) in the .HEAD file of the experiment-aligned Parent Volume */
450 VP->CENTER_BASE = NULL; /*!< pointer to the named attribute (3x1) in the .HEAD file of the experiment-aligned Parent Volume */
451 VP->MATVEC = NULL; /*!< pointer to the named attribute (12x1) in the .HEAD file of the experiment-aligned Parent Volume */
452 VP->Hand = 1; /*!< Handedness of axis 1 RH, -1 LH*/
453 VP->MATVEC_source = NO_WARP;
454 SUMA_RETURN(VP);
455 }
SUMA_Free_VolPar(SUMA_VOLPAR * VP)456 SUMA_Boolean SUMA_Free_VolPar (SUMA_VOLPAR *VP)
457 {
458 static char FuncName[]={"SUMA_Free_VolPar"};
459
460 SUMA_ENTRY;
461 if (!VP) SUMA_RETURN (YUP);
462 if (VP->prefix != NULL) SUMA_free(VP->prefix);
463 if (VP->headname != NULL) SUMA_free(VP->headname);
464 if (VP->idcode_str != NULL) SUMA_free(VP->idcode_str);
465 if (VP->filecode != NULL) SUMA_free(VP->filecode);
466 if (VP->dirname != NULL) SUMA_free(VP->dirname);
467 if (VP->vol_idcode_str != NULL) SUMA_free(VP->vol_idcode_str);
468 if (VP->vol_idcode_date != NULL) SUMA_free(VP->vol_idcode_date);
469 if (VP->CENTER_OLD != NULL) SUMA_free(VP->CENTER_OLD);
470 if (VP->CENTER_BASE != NULL) SUMA_free(VP->CENTER_BASE);
471 if (VP->MATVEC != NULL) SUMA_free(VP->MATVEC);
472 if (VP != NULL) SUMA_free(VP);
473 SUMA_RETURN (YUP);
474 }
475
476
SUMA_VolParFromDset(THD_3dim_dataset * dset)477 SUMA_VOLPAR *SUMA_VolParFromDset (THD_3dim_dataset *dset)
478 {
479 ATR_float *atr=NULL, *atrkeep=NULL;
480 static char FuncName[]={"SUMA_VolParFromDset"};
481 SUMA_VOLPAR *VP=NULL;
482 int ii, nxform = 0;
483 MCW_idcode idcode;
484 SUMA_Boolean LocalHead = NOPE;
485
486 SUMA_ENTRY;
487
488 /* read the header of the parent volume */
489 if (dset == NULL) {
490 fprintf (SUMA_STDERR,"Error %s:\nNULL dset\n", FuncName);
491 SUMA_RETURN (NULL);
492 }
493
494 /* allocate for VP */
495 VP = SUMA_Alloc_VolPar();
496 if (VP == NULL) {
497 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Alloc_VolPar\n", FuncName);
498 SUMA_RETURN (NULL);
499 }
500
501 /* retain pertinent info */
502 VP->isanat = ISANAT(dset);
503 VP->nx = DSET_NX(dset);
504 VP->ny = DSET_NY(dset);
505 VP->nz = DSET_NZ(dset);
506 VP->dx = DSET_DX(dset);
507 VP->dy = DSET_DY(dset);
508 VP->dz = DSET_DZ(dset);
509 VP->xorg = DSET_XORG(dset);
510 VP->yorg = DSET_YORG(dset);
511 VP->zorg = DSET_ZORG(dset);
512 ii = strlen(DSET_HEADNAME(dset));
513 VP->headname = (char *)SUMA_malloc(ii+1);
514 ii = strlen(DSET_PREFIX(dset));
515 VP->prefix = (char *)SUMA_malloc(ii+1);
516 ii = strlen(DSET_FILECODE(dset));
517 VP->filecode = (char *)SUMA_malloc(ii+1);
518 ii = strlen(DSET_DIRNAME(dset));
519 VP->dirname = (char *)SUMA_malloc(ii+1);
520 ii = strlen(dset->idcode.str);
521 VP->vol_idcode_str = (char *)SUMA_malloc(ii+1);
522 ii = strlen(dset->idcode.date);
523 VP->vol_idcode_date = (char *)SUMA_malloc(ii+1);
524 if ( VP->prefix == NULL ||
525 VP->headname == NULL ||
526 VP->filecode == NULL ||
527 VP->vol_idcode_date == NULL ||
528 VP->dirname == NULL ||
529 VP->vol_idcode_str == NULL) {
530 fprintf( SUMA_STDERR,
531 "Error %s: Failed to allocate for strings. Kill me, please.\n",
532 FuncName);
533 SUMA_Free_VolPar(VP);
534 SUMA_RETURN (NULL);
535 }
536 VP->prefix = strcpy(VP->prefix, DSET_PREFIX(dset));
537 VP->headname = strcpy(VP->headname, DSET_HEADNAME(dset));
538 VP->filecode = strcpy(VP->filecode, DSET_FILECODE(dset));
539 VP->dirname = strcpy(VP->dirname, DSET_DIRNAME(dset));
540 VP->vol_idcode_str = strcpy(VP->vol_idcode_str, dset->idcode.str);
541 VP->vol_idcode_date = strcpy(VP->vol_idcode_date, dset->idcode.date);
542 VP->xxorient = dset->daxes->xxorient;
543 VP->yyorient = dset->daxes->yyorient;
544 VP->zzorient = dset->daxes->zzorient;
545
546 if (LocalHead) {
547 fprintf (SUMA_STDERR,"%s: dset->idcode_str = %s\n",
548 FuncName, dset->idcode.str);
549 fprintf (SUMA_STDERR,"%s: VP->vol_idcode_str = %s\n",
550 FuncName, VP->vol_idcode_str);
551 }
552
553 /* Get the 3drotate matrix if possible*/
554 VP->MATVEC_source = NO_WARP;
555 if ((atr = THD_find_float_atr( dset->dblk , "ROTATE_MATVEC_000000" ))) {
556 if (!atrkeep) {
557 atrkeep = atr;
558 VP->MATVEC_source = ROTATE_WARP;
559 }
560 ++nxform;
561 }
562 if ((atr = THD_find_float_atr( dset->dblk , "TAGALIGN_MATVEC" ))) {
563 if (!atrkeep) {
564 atrkeep = atr;
565 VP->MATVEC_source = TAGALIGN_WARP;
566 }
567 ++nxform;
568 }
569 if ((atr = THD_find_float_atr( dset->dblk , "WARPDRIVE_MATVEC_INV_000000" ))){
570 if (!atrkeep) {
571 atrkeep = atr;
572 VP->MATVEC_source = WARPDRIVE_WARP;
573 }
574 ++nxform;
575 }
576 if ((atr = THD_find_float_atr( dset->dblk , "ALLINEATE_MATVEC_S2B_000000" ))){
577 if (!atrkeep) {
578 atrkeep = atr;
579 VP->MATVEC_source = ALLINEATE_WARP;
580 }
581 ++nxform;
582 }
583 if ((atr = THD_find_float_atr( dset->dblk , "VOLREG_MATVEC_000000" ))) {
584 if (!atrkeep) {
585 atrkeep = atr;
586 VP->MATVEC_source = VOLREG_WARP;
587 }
588 ++nxform;
589 }
590 atr = atrkeep;
591 if (nxform > 1) {
592 SUMA_S_Warnv("More than one (%d) plausible transforms found.\n"
593 "Abiding by %s\n",
594 nxform, SUMA_WarpTypeName(VP->MATVEC_source));
595 }
596
597 if (atr == NULL) {
598 VP->MATVEC = NULL;
599 }else {
600 VP->MATVEC = (double *)SUMA_calloc(12, sizeof(double));
601 if (VP->MATVEC != NULL) {
602 if (atr->nfl == 12) {
603 for (ii=0; ii<12; ++ii) VP->MATVEC[ii] = (double)atr->fl[ii];
604 } else {
605 fprintf( SUMA_STDERR,
606 "Error %s: MATVEC does not have 12 elements.\n",
607 FuncName);
608 }
609 } else {
610 fprintf(SUMA_STDERR,
611 "Error %s: Failed to allocate for VP->MATVEC\n",
612 FuncName);
613 }
614 }
615
616
617 /* Get the center base coordinates */
618 switch (VP->MATVEC_source) {
619 case VOLREG_WARP:
620 atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_BASE");
621 break;
622 case ROTATE_WARP:
623 atr = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_BASE");
624 break;
625 default:
626 atr = NULL;
627 break;
628 }
629
630 if (atr == NULL) {
631 VP->CENTER_BASE = NULL;
632 } else {
633 VP->CENTER_BASE = (double *)SUMA_calloc(3, sizeof(double));
634 if (VP->CENTER_BASE != NULL) {
635 if (atr->nfl == 3) {
636 for (ii=0; ii<3; ++ii) VP->CENTER_BASE[ii] = atr->fl[ii];
637 } else {
638 fprintf( SUMA_STDERR,
639 "Error %s: CENTER_BASE does not have 12 elements.\n",
640 FuncName);
641 }
642 } else {
643 fprintf( SUMA_STDERR,
644 "Error %s: Failed to allocate for VP->CENTER_BASE\n",
645 FuncName);
646 }
647 }
648
649 /* CENTER_OLD */
650 switch (VP->MATVEC_source) {
651 case VOLREG_WARP:
652 atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_OLD");
653 break;
654 case ROTATE_WARP:
655 atr = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_OLD");
656 break;
657 default:
658 atr = NULL;
659 break;
660 }
661
662 if (atr == NULL) {
663 VP->CENTER_OLD = NULL;
664 } else {
665 VP->CENTER_OLD = (double *)SUMA_calloc(3, sizeof(double));
666 if (VP->CENTER_OLD != NULL) {
667 if (atr->nfl == 3) {
668 for (ii=0; ii<3; ++ii) VP->CENTER_OLD[ii] = atr->fl[ii];
669 } else {
670 fprintf( SUMA_STDERR,
671 "Error %s: CENTER_OLD does not have 12 elements.\n",
672 FuncName);
673 }
674 } else {
675 fprintf( SUMA_STDERR,
676 "Error %s: Failed to allocate for VP->CENTER_OLD\n", FuncName);
677 }
678 }
679
680 /* handedness */
681 VP->Hand = SUMA_THD_handedness( dset );
682
683 SUMA_RETURN (VP);
684 }
685
SUMA_VolPar_Attr(char * volparent_name)686 SUMA_VOLPAR *SUMA_VolPar_Attr (char *volparent_name)
687 {
688 ATR_float *atr;
689 static char FuncName[]={"SUMA_VolPar_Attr"};
690 SUMA_VOLPAR *VP;
691 THD_3dim_dataset *dset;
692 int ii;
693 MCW_idcode idcode;
694 SUMA_Boolean LocalHead = NOPE;
695
696 SUMA_ENTRY;
697
698 /* read the header of the parent volume */
699 dset = THD_open_dataset(volparent_name);
700 if (dset == NULL) {
701 fprintf (SUMA_STDERR,
702 "Error %s: Could not read %s\n", FuncName, volparent_name);
703 SUMA_RETURN (NULL);
704 }
705
706 VP = SUMA_VolParFromDset(dset);
707 if (!VP) {
708 SUMA_SL_Err("Failed in SUMA_VolParFromDset");
709 }
710
711 /* now free the dset pointer */
712 THD_delete_3dim_dataset( dset , False ) ;
713
714 SUMA_RETURN (VP);
715 }
716
717 /*!
718 \brief Form a string containing the info of the volume parent
719
720 \param VP (SUMA_VOLPAR *) Volume parent structure.
721 \return s (char *) pointer to NULL terminated string containing surface info.
722 It is your responsability to free it.
723
724 \sa SUMA_Show_VolPar
725 */
726
SUMA_VolPar_Info(SUMA_VOLPAR * VP)727 char *SUMA_VolPar_Info (SUMA_VOLPAR *VP)
728 {
729 static char FuncName[]={"SUMA_VolPar_Info"};
730 char stmp[1000], *s=NULL;
731 SUMA_STRING *SS = NULL;
732
733 SUMA_ENTRY;
734
735 SS = SUMA_StringAppend (NULL, NULL);
736
737 if (VP) {
738 sprintf (stmp,"\nVP contents:\n");
739 SS = SUMA_StringAppend (SS, stmp);
740 sprintf (stmp,
741 "prefix: %s\tfilecode: %s\theadname %s\tdirname: %s\n"
742 "Id code str:%s\tID code date: %s\n",
743 VP->prefix, VP->filecode, VP->headname, VP->dirname,
744 VP->vol_idcode_str, VP->vol_idcode_date);
745 SS = SUMA_StringAppend (SS, stmp);
746 if (VP->idcode_str) SS = SUMA_StringAppend (SS, "IDcode is NULL\n");
747 else SS = SUMA_StringAppend_va (SS, "IDcode: %s\n", VP->idcode_str);
748
749 sprintf (stmp,"isanat: %d\n", VP->isanat);
750 SS = SUMA_StringAppend (SS, stmp);
751 sprintf (stmp,"Orientation: %d %d %d\n",
752 VP->xxorient, VP->yyorient, VP->zzorient);
753 if (VP->Hand == 1)
754 SS = SUMA_StringAppend (SS, "Right Hand Coordinate System.\n");
755 else if (VP->Hand == -1)
756 SS = SUMA_StringAppend (SS, "Left Hand Coordinate System.\n");
757 else SS = SUMA_StringAppend (SS, "No hand coordinate system!\n");
758
759 SS = SUMA_StringAppend (SS, stmp);
760 sprintf (stmp,"Origin: %f %f %f\n",
761 VP->xorg, VP->yorg, VP->zorg);
762 SS = SUMA_StringAppend (SS, stmp);
763 sprintf (stmp,"Delta: %f %f %f\n",
764 VP->dx, VP->dy, VP->dz);
765 SS = SUMA_StringAppend (SS, stmp);
766 sprintf (stmp,"N: %d %d %d\n",
767 VP->nx, VP->ny, VP->nz);
768 SS = SUMA_StringAppend (SS, stmp);
769
770 SS = SUMA_StringAppend_va( SS,
771 "VolPar transform type: %d\n",
772 SUMA_WarpTypeName(VP->MATVEC_source));
773 if (VP->MATVEC != NULL) {
774 sprintf (stmp,"VP->MATVEC = \n\tMrot\tDelta\n");
775 SS = SUMA_StringAppend (SS, stmp);
776 sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n",
777 VP->MATVEC[0], VP->MATVEC[1], VP->MATVEC[2], VP->MATVEC[3]);
778 SS = SUMA_StringAppend (SS, stmp);
779 sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n",
780 VP->MATVEC[4], VP->MATVEC[5], VP->MATVEC[6], VP->MATVEC[7]);
781 SS = SUMA_StringAppend (SS, stmp);
782 sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n",
783 VP->MATVEC[8], VP->MATVEC[9], VP->MATVEC[10], VP->MATVEC[11]);
784 SS = SUMA_StringAppend (SS, stmp);
785 } else {
786 sprintf (stmp,"VP->MATVEC = NULL\n");
787 SS = SUMA_StringAppend (SS, stmp);
788 }
789 if (VP->CENTER_OLD != NULL) {
790 sprintf (stmp,"VP->CENTER_OLD = %f, %f, %f\n",
791 VP->CENTER_OLD[0], VP->CENTER_OLD[1], VP->CENTER_OLD[2]);
792 SS = SUMA_StringAppend (SS, stmp);
793 }else {
794 sprintf (stmp,"VP->CENTER_OLD = NULL\n");
795 SS = SUMA_StringAppend (SS, stmp);
796 }
797 if (VP->CENTER_BASE != NULL) {
798 sprintf (stmp,"VP->CENTER_BASE = %f, %f, %f\n",
799 VP->CENTER_BASE[0], VP->CENTER_BASE[1], VP->CENTER_BASE[2]);
800 SS = SUMA_StringAppend (SS, stmp);
801 } else {
802 sprintf (stmp,"VP->CENTER_BASE = NULL\n");
803 SS = SUMA_StringAppend (SS, stmp);
804 }
805 }else{
806 sprintf (stmp, "NULL Volume Parent Pointer.\n");
807 SS = SUMA_StringAppend (SS, stmp);
808 }
809
810 /* clean SS */
811 SS = SUMA_StringAppend (SS, NULL);
812 /* copy s pointer and free SS */
813 s = SS->s;
814 SUMA_free(SS);
815
816 SUMA_RETURN (s);
817 }
818 /*!
819 Show the contents of SUMA_VOLPAR structure
820 \sa SUMA_VolPar_Info
821 */
SUMA_Show_VolPar(SUMA_VOLPAR * VP,FILE * Out)822 void SUMA_Show_VolPar(SUMA_VOLPAR *VP, FILE *Out)
823 {
824 static char FuncName[]={"SUMA_Show_VolPar"};
825 char *s;
826
827 SUMA_ENTRY;
828
829 if (Out == NULL) Out = SUMA_STDOUT;
830
831 s = SUMA_VolPar_Info(VP);
832
833 if (s) {
834 fprintf (Out, "%s", s);
835 SUMA_free(s);
836 }else {
837 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_VolPar_Info.\n", FuncName);
838 }
839
840 SUMA_RETURNe;
841
842 }
843
844 /*!
845 \brief Transform the coordinates of a surface object to AFNI-DICOM convention and transform the coordinates by SO->VolPar
846 ans = SUMA_Align_to_VolPar (SO, SF_Struct);
847 \param SO (SUMA_SurfaceObject *)
848 \param S_Struct (void *) That is only needed for SureFit surfaces and is nothing but a type cast of a SUMA_SureFit_struct containing information on cropping.
849 send NULL for all other surfaces.
850 \return YUP/NOPE
851 For SureFit and FreeSurfer surfaces, the coordinates are first set in RAI (DICOM) coordinate system before applying SO->VolPar.
852 For other surface formats, SO->VolPar is applied to whatever coordinates are in SO->NodeList
853 */
SUMA_Align_to_VolPar(SUMA_SurfaceObject * SO,void * S_Struct)854 SUMA_Boolean SUMA_Align_to_VolPar (SUMA_SurfaceObject *SO, void * S_Struct)
855 {
856 static char FuncName[]={"SUMA_Align_to_VolPar"};
857 char orcode[6];
858 float xx, yy, zz;
859 THD_coorder * cord_surf, *cord_RAI;
860 int i, ND, id;
861 SUMA_SureFit_struct *SF;
862 SUMA_FreeSurfer_struct *FS;
863 SUMA_Boolean LocalHead = NOPE;
864
865 SUMA_ENTRY;
866
867 SUMA_LHv("Working %s\n",SO->Name.FileName);
868
869 /* don't bother if surface is flat */
870 if (SUMA_is_Constant_Z_Coord(SO->NodeList, SO->N_Node, 0.001)) {
871 /* you can also use the more generic :
872 SUMA_is_Flat_Surf_Coords_PCA(SO->NodeList, SO->N_Node, 0.001, 0.01)
873 but that one's slower */
874 SUMA_LH("This one is flat\n");
875 SO->APPLIED_A2Exp_XFORM = NO_WARP;
876 SUMA_RETURN(YUP);
877 }
878 if (SO->isSphere == SUMA_GEOM_NOT_SET) SUMA_SetSphereParams(SO, -0.1);
879 if (SUMA_IS_GEOM_SYMM(SO->isSphere)) { /* March 2013 */
880 SUMA_LH("This one is a ball\n");
881 SO->APPLIED_A2Exp_XFORM = NO_WARP;
882 SUMA_RETURN(YUP);
883 }
884 /* allocate for cord structure */
885 cord_surf = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
886 cord_RAI = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
887 if (cord_surf == NULL || cord_RAI == NULL) {
888 fprintf(SUMA_STDERR,"Error %s: failed to allocate for cord\n", FuncName);
889 SUMA_RETURN (NOPE);
890 }
891
892 /* do the dance to get cord for RAI */
893 THD_coorder_fill( NULL, cord_RAI);
894 ND = SO->NodeDim;
895 switch (SO->FileType) {
896 case SUMA_INVENTOR_GENERIC:
897 case SUMA_OPENDX_MESH:
898 case SUMA_OBJ_MESH:
899 case SUMA_PREDEFINED:
900 case SUMA_PLY:
901 case SUMA_STL:
902 case SUMA_BYU:
903 case SUMA_VEC:
904 /* Do nothing */
905 break;
906 case SUMA_FREE_SURFER:
907 case SUMA_FREE_SURFER_PATCH:
908 /* For free surfer, all you need to do is
909 go from LPI to RAI (DICOM)*/
910 sprintf(orcode,"LPI");
911 THD_coorder_fill(orcode , cord_surf);
912 /*loop over XYZs and change them to dicom*/
913 for (i=0; i < SO->N_Node; ++i) {
914 id = i * ND;
915 THD_coorder_to_dicom (cord_surf,
916 &(SO->NodeList[id]),
917 &(SO->NodeList[id+1]),
918 &(SO->NodeList[id+2]));
919 }
920 break;
921 case SUMA_SUREFIT:
922 /* For SureFit, coordinates are actually a
923 float version of the indices */
924 SF = (SUMA_SureFit_struct *)S_Struct;
925 if (SF->caret_version < 5.2) {
926 THD_fvec3 fv, iv;
927 float D[3];
928 /* Calcluate Delta caused by cropping */
929 for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
930 SUMA_LHv("caret_version: %f\n"
931 "AC_WholeVolume: [%f %f %f]\n"
932 "AC: [%f %f %f]\n"
933 "Shift Values: [%f, %f, %f]\n"
934 "Node 0 init: [%f, %f, %f]\n"
935 "CropMin: [%f, %f, %f]\n"
936 "CropMax: [%f, %f, %f]\n",
937 SF->caret_version,
938 SF->AC_WholeVolume[0],SF->AC_WholeVolume[1],
939 SF->AC_WholeVolume[2],
940 SF->AC[0], SF->AC[1], SF->AC[2],
941 D[0], D[1], D[2],
942 SO->NodeList[0], SO->NodeList[1],SO->NodeList[2],
943 SF->CropMin[0], SF->CropMin[1], SF->CropMin[2],
944 SF->CropMax[0], SF->CropMax[1], SF->CropMax[2]);
945 for (i=0; i < SO->N_Node; ++i) {
946 id = i * ND;
947 /* change float indices to mm coords */
948 iv.xyz[0] = SO->NodeList[id] + D[0];
949 iv.xyz[1] = SO->NodeList[id+1] + D[1];
950 iv.xyz[2] = SO->NodeList[id+2] + D[2];
951 fv = SUMA_THD_3dfind_to_3dmm( SO, iv );
952
953 /* change mm to RAI coords */
954 iv = SUMA_THD_3dmm_to_dicomm( SO->VolPar->xxorient,
955 SO->VolPar->yyorient,
956 SO->VolPar->zzorient, fv );
957 SO->NodeList[id] = iv.xyz[0];
958 SO->NodeList[id+1] = iv.xyz[1];
959 SO->NodeList[id+2] = iv.xyz[2];
960 }
961 SUMA_LHv("Node 0 RAI: [%f, %f, %f]\n",
962 SO->NodeList[0], SO->NodeList[1],SO->NodeList[2]);
963 } else {
964 float D[3]={0.0, 0.0, 0.0};
965 /* Calcluate Delta caused by cropping */
966 for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
967 SUMA_LHv("caret_version: %f\n"
968 "AC_WholeVolume: [%f %f %f]\n"
969 "AC: [%f %f %f]\n"
970 "Shift Values: [%f, %f, %f]\n"
971 "Node 0 init: [%f, %f, %f]\n"
972 "CropMin: [%f, %f, %f]\n"
973 "CropMax: [%f, %f, %f]\n",
974 SF->caret_version,
975 SF->AC_WholeVolume[0],SF->AC_WholeVolume[1],
976 SF->AC_WholeVolume[2],
977 SF->AC[0], SF->AC[1], SF->AC[2],
978 D[0], D[1], D[2],
979 SO->NodeList[0], SO->NodeList[1],SO->NodeList[2],
980 SF->CropMin[0], SF->CropMin[1], SF->CropMin[2],
981 SF->CropMax[0], SF->CropMax[1], SF->CropMax[2]);
982 if ( SUMA_ABS(D[0]) > 0.0000001 ||
983 SUMA_ABS(D[1]) > 0.0000001 ||
984 SUMA_ABS(D[2]) > 0.0000001 ) {
985 SUMA_S_Notev("Shift Values: [%f, %f, %f]\n",
986 D[0], D[1], D[2]);
987 SUMA_S_Note( "If surface alignment is off,\n"
988 "notify authors and send sample data.\n");
989
990 }
991 /* Caret, just LPI baby, take it to RAI*/
992 for (i=0; i < SO->N_Node; ++i) {
993 id = i * ND;
994 SO->NodeList[id] = -SO->NodeList[id];
995 SO->NodeList[id+1] = -SO->NodeList[id+1];
996 }
997
998 SUMA_LHv("Node 0 RAI : [%f, %f, %f]\n",
999 SO->NodeList[0], SO->NodeList[1],SO->NodeList[2]);
1000 }
1001 break;
1002 case SUMA_BRAIN_VOYAGER:
1003 #if 0
1004 /* this contraption was used when BV's afni volumes had 0,0,0 for origin
1005 which is inappropriate */
1006 /* For Brain Voyager, all you need to do is
1007 go from ASR to RAI (DICOM)
1008 Note: The center of the volume is at the 1st voxel's center and
1009 that huge
1010 center shift, relative to standard AFNI dsets (centered about
1011 middle of volume)
1012 might throw off 3dVolreg. If you want to shift volume's center to be in
1013 the middle voxel, you'll need to shift the surface coordinates
1014 before transforming
1015 them to RAI*/
1016 sprintf(orcode,"ASR");
1017 THD_coorder_fill(orcode , cord_surf);
1018 /*loop over XYZs and change them to dicom*/
1019 for (i=0; i < SO->N_Node; ++i) {
1020 id = i * ND;
1021 THD_coorder_to_dicom (cord_surf,
1022 &(SO->NodeList[id]),
1023 &(SO->NodeList[id+1]),
1024 &(SO->NodeList[id+2]));
1025 }
1026 #else /*ZSS: Nov. 1 07 */
1027 if (SO->VolPar) {
1028 /* looks like coordinates are in float index units, go to dicomm */
1029 if (!SUMA_vec_3dfind_to_dicomm(SO->NodeList,SO->N_Node,SO->VolPar)) {
1030 SUMA_S_Err("Failed to xform coords.");
1031 SUMA_RETURN (NOPE);
1032 }
1033 }
1034 #endif
1035 break;
1036 case SUMA_GIFTI: /* have to apply coord xform */
1037
1038 break;
1039 default:
1040 fprintf( SUMA_STDERR,
1041 "Warning %s: Unknown SO->FileType.\n"
1042 "Assuming coordinates are in DICOM already.\n", FuncName);
1043 break;
1044 }
1045
1046 if (SO->VolPar && !SUMA_Apply_VolReg_Trans (SO)) {
1047 fprintf( SUMA_STDERR,
1048 "Error %s: Failed in SUMA_Apply_VolReg_Trans.\n", FuncName);
1049 SUMA_RETURN (NOPE);
1050 }
1051
1052 SUMA_nixSODim(SO); /* Dims need recomputing, this might be overkill if
1053 nothing was done above */
1054 SUMA_SetSODims(SO);
1055 SUMA_RETURN (YUP);
1056 }
1057
1058 /*!
1059 \brief Undo the transform the coordinates of a surface object from AFNI-DICOM convention to the native space
1060 Transforms in SO->VolPar are ignored
1061 ans = SUMA_Delign_to_VolPar (SO, SF_Struct);
1062 \param SO (SUMA_SurfaceObject *)
1063 \param S_Struct (void *) That is only needed for SureFit surfaces and is nothing but a type cast of a SUMA_SureFit_struct containing information on cropping.
1064 send NULL for all other surfaces.
1065 \return YUP/NOPE
1066 For SureFit and FreeSurfer surfaces, the coordinates are first set in RAI (DICOM) coordinate system before applying SO->VolPar.
1067 For other surface formats, SO->VolPar is applied to whatever coordinates are in SO->NodeList
1068 */
SUMA_Delign_to_VolPar(SUMA_SurfaceObject * SO,void * S_Struct)1069 SUMA_Boolean SUMA_Delign_to_VolPar (SUMA_SurfaceObject *SO, void * S_Struct)
1070 {
1071 static char FuncName[]={"SUMA_Delign_to_VolPar"};
1072 char orcode[6];
1073 float xx, yy, zz;
1074 THD_coorder * cord_surf, *cord_RAI;
1075 int i, ND, id;
1076 SUMA_SureFit_struct *SF;
1077 SUMA_FreeSurfer_struct *FS;
1078
1079 SUMA_ENTRY;
1080
1081 /* allocate for cord structure */
1082 cord_surf = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
1083 cord_RAI = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
1084 if (cord_surf == NULL || cord_RAI == NULL) {
1085 fprintf(SUMA_STDERR,"Error %s: failed to allocate for cord\n", FuncName);
1086 SUMA_RETURN (NOPE);
1087 }
1088
1089 /* do the dance to get cord for RAI */
1090 THD_coorder_fill( NULL, cord_RAI);
1091 ND = SO->NodeDim;
1092 switch (SO->FileType) {
1093 case SUMA_INVENTOR_GENERIC:
1094 case SUMA_OPENDX_MESH:
1095 case SUMA_OBJ_MESH:
1096 case SUMA_PREDEFINED:
1097 case SUMA_PLY:
1098 case SUMA_STL:
1099 case SUMA_BYU:
1100 case SUMA_VEC:
1101 /* Do nothing */
1102 break;
1103 case SUMA_FREE_SURFER:
1104 case SUMA_FREE_SURFER_PATCH:
1105 /* For free surfer, all you need to do is
1106 go from LPI to RAI (DICOM)*/
1107 sprintf(orcode,"LPI");
1108 THD_coorder_fill(orcode , cord_surf);
1109 /*loop over XYZs and change them to surface's order*/
1110 for (i=0; i < SO->N_Node; ++i) {
1111 id = i * ND;
1112 THD_dicom_to_coorder (cord_surf, &(SO->NodeList[id]), &(SO->NodeList[id+1]), &(SO->NodeList[id+2]));
1113 }
1114 break;
1115 case SUMA_SUREFIT:
1116 SF = (SUMA_SureFit_struct *)S_Struct;
1117 if (SF->caret_version < 5.2) {
1118 /* For SureFit, coordinates are actually a float version of the indices */
1119 SUMA_SL_Warn( "Reverse implementation not finished\n"
1120 "Send me a complaint, I must have forgotten\n"
1121 "Coords will be left untouched. (saadz@mail.nih.gov)\n");
1122 #if 0
1123 { THD_fvec3 fv, iv;
1124 float D[3];
1125 /* Calcluate Delta caused by cropping */
1126 for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
1127 /* fprintf (SUMA_STDERR,"%s: Shift Values: [%f, %f, %f]\n", FuncName, D[0], D[1], D[2]); */
1128 for (i=0; i < SO->N_Node; ++i) {
1129 id = i * ND;
1130 /* change float indices to mm coords */
1131 iv.xyz[0] = SO->NodeList[id] + D[0];
1132 iv.xyz[1] = SO->NodeList[id+1] + D[1];
1133 iv.xyz[2] = SO->NodeList[id+2] + D[2];
1134 fv = SUMA_THD_3dfind_to_3dmm( SO, iv );
1135
1136 /* change mm to RAI coords */
1137 iv = SUMA_THD_3dmm_to_dicomm( SO->VolPar->xxorient, SO->VolPar->yyorient, SO->VolPar->zzorient, fv );
1138 SO->NodeList[id] = iv.xyz[0];
1139 SO->NodeList[id+1] = iv.xyz[1];
1140 SO->NodeList[id+2] = iv.xyz[2];
1141 }
1142 }
1143 #endif
1144 } else {
1145 float D[3];
1146 /* Calcluate Delta caused by cropping */
1147 for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
1148 if (D[0] != 0.0f || D[1] != 0.0f ||D[2] != 0.0f) {
1149 fprintf (SUMA_STDERR,"Error %s: Shift Values: [%f, %f, %f]\n", FuncName, D[0], D[1], D[2]);
1150 fprintf (SUMA_STDERR,"Never encountered this case. Please notify authors and send sample data.\n");
1151 SUMA_RETURN (NOPE);
1152 }
1153
1154 /* just go back from RAI to LPI */
1155 for (i=0; i < SO->N_Node; ++i) {
1156 id = i * ND;
1157 SO->NodeList[id] = -SO->NodeList[id];
1158 SO->NodeList[id+1] = -SO->NodeList[id+1];
1159 }
1160 }
1161 break;
1162 case SUMA_BRAIN_VOYAGER:
1163 #if 0
1164 /* For Brain Voyager, all you need to do is
1165 go from AIR to RAI (DICOM)
1166 Note: The center of the volume is at the 1st voxel's center and that huge
1167 center shift, relative to standard AFNI dsets (centered about middle of volume)
1168 might throw off 3dVolreg. If you want to shift volume's center to be in
1169 the middle voxel, you'll need to shift the surface coordinates before transforming
1170 them to ASR*/
1171 sprintf(orcode,"ASR");
1172 THD_coorder_fill(orcode , cord_surf);
1173 /*loop over XYZs and change them to native space*/
1174 for (i=0; i < SO->N_Node; ++i) {
1175 id = i * ND;
1176 THD_dicom_to_coorder (cord_surf, &(SO->NodeList[id]), &(SO->NodeList[id+1]), &(SO->NodeList[id+2]));
1177 }
1178 #else /*ZSS: Nov. 1 07 */
1179 if (SO->VolPar) {
1180 /* to back to float index units */
1181 if (!SUMA_vec_dicomm_to_3dfind(SO->NodeList, SO->N_Node, SO->VolPar)) {
1182 SUMA_S_Err("Failed to xform coords.");
1183 SUMA_RETURN (NOPE);
1184 }
1185 } else {
1186 fprintf(SUMA_STDERR,"Error %s: Can't delign witout a volpar.\n", FuncName);
1187 SUMA_RETURN (NOPE);
1188 }
1189 #endif
1190 break;
1191 default:
1192 fprintf(SUMA_STDERR,"Warning %s: Unknown SO->FileType. Assuming coordinates are in DICOM already.\n", FuncName);
1193 break;
1194 }
1195
1196 #if 0
1197 /* I don't think the inverse of that step will be needed .... */
1198 if (!SUMA_Apply_VolReg_Trans (SO)) {
1199 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Apply_VolReg_Trans.\n", FuncName);
1200 SUMA_RETURN (NOPE);
1201 }
1202 #endif
1203
1204 SUMA_nixSODim(SO); /* Dims need recomputing, this might be overkill if
1205 nothing was done above */
1206 SUMA_SetSODims(SO);
1207 SUMA_RETURN (YUP);
1208 }
1209
1210
SUMA_Apply_Coord_xform(float * NodeList,int N_Node,int NodeDim,double Xform[4][4],int doinv,double * ppshift)1211 SUMA_Boolean SUMA_Apply_Coord_xform(float *NodeList,
1212 int N_Node,
1213 int NodeDim,
1214 double Xform[4][4],
1215 int doinv,
1216 double *ppshift)
1217 {
1218 static char FuncName[]={"SUMA_Apply_Coord_xform"};
1219 double x, y, z;
1220 int i=0, id = 0;
1221 mat44 A, A0;
1222 SUMA_Boolean LocalHead = NOPE;
1223
1224 SUMA_ENTRY;
1225
1226 if (!NodeList) SUMA_RETURN(NOPE);
1227
1228 /* check for identity */
1229 if ( SUMA_IS_XFORM_IDENTITY(Xform) ) {
1230 SUMA_LH("Indentity, nothing to do.");
1231 SUMA_RETURN(YUP);
1232 }
1233
1234 if (!doinv) {
1235 LOAD_MAT44( A, \
1236 Xform[0][0], Xform[0][1], Xform[0][2], Xform[0][3], \
1237 Xform[1][0], Xform[1][1], Xform[1][2], Xform[1][3], \
1238 Xform[2][0], Xform[2][1], Xform[2][2], Xform[2][3] );
1239 } else {
1240 SUMA_LH("Inverting xform");
1241 LOAD_MAT44( A0, \
1242 Xform[0][0], Xform[0][1], Xform[0][2], Xform[0][3], \
1243 Xform[1][0], Xform[1][1], Xform[1][2], Xform[1][3], \
1244 Xform[2][0], Xform[2][1], Xform[2][2], Xform[2][3] );
1245 A = nifti_mat44_inverse(A0);
1246 }
1247
1248 SUMA_LHv("Node 0, Pre: %f %f %f\n",
1249 NodeList[0] , NodeList[1], NodeList[2]);
1250 for (i=0; i < N_Node; ++i) {
1251 id = NodeDim * i;
1252 x = (double)NodeList[id] ;
1253 y = (double)NodeList[id+1] ;
1254 z = (double)NodeList[id+2] ;
1255 if (ppshift) {
1256 x += ppshift[0];
1257 y += ppshift[1];
1258 z += ppshift[2];
1259 }
1260
1261 /* Apply the rotation matrix XYZn = Mrot x XYZ + Delta*/
1262 NodeList[id] = (float) ( A.m[0][0] * x +
1263 A.m[0][1] * y +
1264 A.m[0][2] * z +
1265 A.m[0][3] );
1266 NodeList[id+1] = (float) ( A.m[1][0] * x +
1267 A.m[1][1] * y +
1268 A.m[1][2] * z +
1269 A.m[1][3] );
1270 NodeList[id+2] = (float) ( A.m[2][0] * x +
1271 A.m[2][1] * y +
1272 A.m[2][2] * z +
1273 A.m[2][3] );
1274 if (ppshift) {
1275 NodeList[id ] -= ppshift[0];
1276 NodeList[id+1] -= ppshift[1];
1277 NodeList[id+2] -= ppshift[2];
1278 }
1279 }
1280 SUMA_LHv("Node 0, Post: %f %f %f\n",
1281 NodeList[0] , NodeList[1], NodeList[2]);
1282
1283 SUMA_RETURN(YUP);
1284 }
1285
1286 /*!
1287 \brief applies the transform in SO->VolPar to SO->NodeList. This only makes sense if
1288 SO->NodeList is in DICOM to begin with...
1289 \param SO (SUMA_SurfaceObject *) You better have NodeList and VolPar in there...
1290 \return YUP/NOPE
1291 \sa SUMA_Align_to_VolPar
1292 NOTE: The field SO->VOLREG_APPLIED is set to NOPE if this function fails, YUP otherwise
1293 */
1294
SUMA_Apply_VolReg_Trans(SUMA_SurfaceObject * SO)1295 SUMA_Boolean SUMA_Apply_VolReg_Trans (SUMA_SurfaceObject *SO)
1296 {
1297 static char FuncName[]={"SUMA_Apply_VolReg_Trans"};
1298 int i, ND, id;
1299 double xform[3][4];
1300 SUMA_Boolean Bad=YUP;
1301 SUMA_Boolean LocalHead=NOPE;
1302
1303 SUMA_ENTRY;
1304
1305 if (get_IgnoreXforms()) {
1306 SUMA_SL_Note("Ignoring any Volreg, TagAlign, Rotate, \n"
1307 "WarpDrive, or Allineate transforms present\n"
1308 "in Surface Volume.\n");
1309 SUMAg_CF->IgnoreVolreg = YUP;
1310 SO->APPLIED_A2Exp_XFORM = NO_WARP;
1311 SUMA_RETURN (YUP);
1312 }
1313
1314 if (!SO || !SO->VolPar) {
1315 SUMA_S_Errv("Bad params, SO=%p, SO->VolPar=%p\n",
1316 SO, SO ? SO->VolPar:NULL);
1317 SUMA_RETURN (NOPE);
1318 }
1319
1320 if (SO->APPLIED_A2Exp_XFORM != 0) {
1321 fprintf (SUMA_STDERR,
1322 "Error %s: \n"
1323 "Volreg (or Tagalign or rotate or warpdrive or allineate)\n"
1324 "already applied. Nothing done.\n", FuncName);
1325 SUMA_RETURN (NOPE);
1326 }
1327
1328 ND = SO->NodeDim;
1329
1330 /* perform the rotation needed to align
1331 the surface with the current experiment's data */
1332 SUMA_LHv("MATVEC_source = %s\n",
1333 SUMA_WarpTypeName(SO->VolPar->MATVEC_source));
1334 switch (SO->VolPar->MATVEC_source) {
1335 case ROTATE_WARP:
1336 case VOLREG_WARP:
1337 if ( SO->VolPar->MATVEC != NULL ||
1338 SO->VolPar->CENTER_OLD != NULL ||
1339 SO->VolPar->CENTER_BASE != NULL) {
1340 Bad = NOPE;
1341 if (SO->VolPar->MATVEC == NULL) {
1342 fprintf(SUMA_STDERR,
1343 "Error %s: SO->VolPar->VOLREG_MATVEC = NULL. \n"
1344 "Cannot perform alignment.\n", FuncName);
1345 Bad = YUP;
1346 }
1347 if (SO->VolPar->CENTER_OLD == NULL) {
1348 fprintf(SUMA_STDERR,
1349 "Error %s: SO->VolPar->CENTER_OLD = NULL.\n"
1350 "Cannot perform alignment.\n", FuncName);
1351 Bad = YUP;
1352 }
1353 if (SO->VolPar->CENTER_BASE == NULL) {
1354 fprintf( SUMA_STDERR,
1355 "Error %s: SO->VolPar->CENTER_BASE = NULL. \n"
1356 "Cannot perform alignment.\n", FuncName);
1357 Bad = YUP;
1358 }
1359 }
1360 break;
1361 case TAGALIGN_WARP:
1362 case ALLINEATE_WARP:
1363 case WARPDRIVE_WARP:
1364 /* all good */
1365 if ( SO->VolPar->MATVEC != NULL) {
1366 Bad = NOPE;
1367 } else {
1368 fprintf(SUMA_STDERR,
1369 "Error %s: SO->VolPar->MATVEC = NULL. \n"
1370 "Cannot perform alignment.\n", FuncName);
1371 Bad = YUP;
1372 }
1373 break;
1374 case NO_WARP:
1375 SO->APPLIED_A2Exp_XFORM = NO_WARP;
1376 SUMA_RETURN (YUP);
1377 default:
1378 fprintf( SUMA_STDERR,
1379 "Error %s: Warptype %d unaccounted for.\n",
1380 FuncName, SO->VolPar->MATVEC_source);
1381 Bad = YUP;
1382 break;
1383 }
1384
1385 /* Now do the transformation */
1386
1387 #if 1
1388 if ( SO->VolPar->MATVEC_source == VOLREG_WARP ||
1389 SO->VolPar->MATVEC_source == ROTATE_WARP ) {
1390 /* remove old center */
1391 for (i=0; i < SO->N_Node; ++i) {
1392 id = ND * i;
1393 SO->NodeList[id ] -= SO->VolPar->CENTER_OLD[0];
1394 SO->NodeList[id+1] -= SO->VolPar->CENTER_OLD[1];
1395 SO->NodeList[id+2] -= SO->VolPar->CENTER_OLD[2];
1396 }
1397 }
1398
1399 /* Now apply affine */
1400 SUMA_Xform1Dto2D(SO->VolPar->MATVEC, xform);
1401 if (!SUMA_Apply_Coord_xform(SO->NodeList, SO->N_Node, SO->NodeDim,
1402 xform, 0, NULL)) {
1403 SUMA_S_Err("Failed to apply AlndExp transform");
1404 Bad = YUP;
1405 }
1406
1407 if ( SO->VolPar->MATVEC_source == VOLREG_WARP ||
1408 SO->VolPar->MATVEC_source == ROTATE_WARP ) {
1409 /* put back new center */
1410 for (i=0; i < SO->N_Node; ++i) {
1411 id = ND * i;
1412 SO->NodeList[id ] += SO->VolPar->CENTER_BASE[0];
1413 SO->NodeList[id+1] += SO->VolPar->CENTER_BASE[1];
1414 SO->NodeList[id+2] += SO->VolPar->CENTER_BASE[2];
1415 }
1416 }
1417
1418 SO->APPLIED_A2Exp_XFORM = SO->VolPar->MATVEC_source;
1419 Bad = NOPE;
1420 #else
1421 /* the olde way, burn it a while after March 27 2008*/
1422 if (UseTagAlign) {
1423 float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
1424
1425 /* fillerup*/
1426 Mrot[0][0] = SO->VolPar->TAGALIGN_MATVEC[0];
1427 Mrot[0][1] = SO->VolPar->TAGALIGN_MATVEC[1];
1428 Mrot[0][2] = SO->VolPar->TAGALIGN_MATVEC[2];
1429 Delta[0] = SO->VolPar->TAGALIGN_MATVEC[3];
1430 Mrot[1][0] = SO->VolPar->TAGALIGN_MATVEC[4];
1431 Mrot[1][1] = SO->VolPar->TAGALIGN_MATVEC[5];
1432 Mrot[1][2] = SO->VolPar->TAGALIGN_MATVEC[6];
1433 Delta[1] = SO->VolPar->TAGALIGN_MATVEC[7];
1434 Mrot[2][0] = SO->VolPar->TAGALIGN_MATVEC[8];
1435 Mrot[2][1] = SO->VolPar->TAGALIGN_MATVEC[9];
1436 Mrot[2][2] = SO->VolPar->TAGALIGN_MATVEC[10];
1437 Delta[2] = SO->VolPar->TAGALIGN_MATVEC[11];
1438
1439 NetShift[0] = Delta[0];
1440 NetShift[1] = Delta[1];
1441 NetShift[2] = Delta[2];
1442
1443 /*
1444 fprintf (SUMA_STDERR,"%s: Applying xform.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
1445 Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
1446 Delta[0], Delta[1], Delta[2]);
1447
1448 */
1449
1450 for (i=0; i < SO->N_Node; ++i) {
1451 id = ND * i;
1452 /* zero the center */
1453 x = SO->NodeList[id] ;
1454 y = SO->NodeList[id+1] ;
1455 z = SO->NodeList[id+2] ;
1456
1457 /* Apply the rotation matrix XYZn = Mrot x XYZ*/
1458 SO->NodeList[id] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
1459 SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
1460 SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
1461
1462 /*apply netshift*/
1463 SO->NodeList[id] += NetShift[0];
1464 SO->NodeList[id+1] += NetShift[1];
1465 SO->NodeList[id+2] += NetShift[2];
1466 }
1467 SO->TAGALIGN_APPLIED = YUP;
1468 } else
1469 SO->TAGALIGN_APPLIED = NOPE;
1470
1471 if (UseWarpAlign) {
1472 float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
1473
1474 /* fillerup*/
1475 Mrot[0][0] = SO->VolPar->WARPDRIVE_MATVEC[0];
1476 Mrot[0][1] = SO->VolPar->WARPDRIVE_MATVEC[1];
1477 Mrot[0][2] = SO->VolPar->WARPDRIVE_MATVEC[2];
1478 Delta[0] = SO->VolPar->WARPDRIVE_MATVEC[3];
1479 Mrot[1][0] = SO->VolPar->WARPDRIVE_MATVEC[4];
1480 Mrot[1][1] = SO->VolPar->WARPDRIVE_MATVEC[5];
1481 Mrot[1][2] = SO->VolPar->WARPDRIVE_MATVEC[6];
1482 Delta[1] = SO->VolPar->WARPDRIVE_MATVEC[7];
1483 Mrot[2][0] = SO->VolPar->WARPDRIVE_MATVEC[8];
1484 Mrot[2][1] = SO->VolPar->WARPDRIVE_MATVEC[9];
1485 Mrot[2][2] = SO->VolPar->WARPDRIVE_MATVEC[10];
1486 Delta[2] = SO->VolPar->WARPDRIVE_MATVEC[11];
1487
1488 NetShift[0] = Delta[0];
1489 NetShift[1] = Delta[1];
1490 NetShift[2] = Delta[2];
1491
1492
1493 /* fprintf (SUMA_STDERR,"%s: Applying xform.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
1494 Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
1495 Delta[0], Delta[1], Delta[2]); */
1496
1497
1498
1499 for (i=0; i < SO->N_Node; ++i) {
1500 id = ND * i;
1501 /* zero the center */
1502 x = SO->NodeList[id] ;
1503 y = SO->NodeList[id+1] ;
1504 z = SO->NodeList[id+2] ;
1505
1506 /* Apply the rotation matrix XYZn = Mrot x XYZ*/
1507 SO->NodeList[id] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
1508 SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
1509 SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
1510
1511 /*apply netshift*/
1512 SO->NodeList[id] += NetShift[0];
1513 SO->NodeList[id+1] += NetShift[1];
1514 SO->NodeList[id+2] += NetShift[2];
1515 }
1516 SO->WARPDRIVE_APPLIED = YUP;
1517 } else
1518 SO->WARPDRIVE_APPLIED = NOPE;
1519
1520 if (UseVolreg) {
1521 float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
1522
1523 /* fillerup*/
1524 Mrot[0][0] = SO->VolPar->VOLREG_MATVEC[0];
1525 Mrot[0][1] = SO->VolPar->VOLREG_MATVEC[1];
1526 Mrot[0][2] = SO->VolPar->VOLREG_MATVEC[2];
1527 Delta[0] = SO->VolPar->VOLREG_MATVEC[3];
1528 Mrot[1][0] = SO->VolPar->VOLREG_MATVEC[4];
1529 Mrot[1][1] = SO->VolPar->VOLREG_MATVEC[5];
1530 Mrot[1][2] = SO->VolPar->VOLREG_MATVEC[6];
1531 Delta[1] = SO->VolPar->VOLREG_MATVEC[7];
1532 Mrot[2][0] = SO->VolPar->VOLREG_MATVEC[8];
1533 Mrot[2][1] = SO->VolPar->VOLREG_MATVEC[9];
1534 Mrot[2][2] = SO->VolPar->VOLREG_MATVEC[10];
1535 Delta[2] = SO->VolPar->VOLREG_MATVEC[11];
1536
1537 NetShift[0] = SO->VolPar->VOLREG_CENTER_BASE[0] + Delta[0];
1538 NetShift[1] = SO->VolPar->VOLREG_CENTER_BASE[1] + Delta[1];
1539 NetShift[2] = SO->VolPar->VOLREG_CENTER_BASE[2] + Delta[2];
1540
1541 /*
1542 fprintf (SUMA_STDERR,"%s: Applying Rotation.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
1543 Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
1544 Delta[0], Delta[1], Delta[2]);
1545 fprintf (SUMA_STDERR,"VOLREG_CENTER_BASE = [%f %f %f]. VOLREG_CENTER_OLD = [%f %f %f]\n", \
1546 SO->VolPar->VOLREG_CENTER_BASE[0], SO->VolPar->VOLREG_CENTER_BASE[1], SO->VolPar->VOLREG_CENTER_BASE[2], \
1547 SO->VolPar->VOLREG_CENTER_OLD[0], SO->VolPar->VOLREG_CENTER_OLD[1], SO->VolPar->VOLREG_CENTER_OLD[2]);
1548 */
1549
1550 for (i=0; i < SO->N_Node; ++i) {
1551 id = ND * i;
1552 /* zero the center */
1553 x = SO->NodeList[id ] - SO->VolPar->VOLREG_CENTER_OLD[0];
1554 y = SO->NodeList[id+1] - SO->VolPar->VOLREG_CENTER_OLD[1];
1555 z = SO->NodeList[id+2] - SO->VolPar->VOLREG_CENTER_OLD[2];
1556
1557 /* Apply the rotation matrix XYZn = Mrot x XYZ*/
1558 SO->NodeList[id ] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
1559 SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
1560 SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
1561
1562 /*apply netshift*/
1563 SO->NodeList[id ] += NetShift[0];
1564 SO->NodeList[id+1] += NetShift[1];
1565 SO->NodeList[id+2] += NetShift[2];
1566 }
1567 SO->VOLREG_APPLIED = YUP;
1568 } else
1569 SO->VOLREG_APPLIED = NOPE;
1570
1571 if (UseRotate) {
1572 float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
1573
1574 /* fillerup*/
1575 Mrot[0][0] = SO->VolPar->ROTATE_MATVEC[0];
1576 Mrot[0][1] = SO->VolPar->ROTATE_MATVEC[1];
1577 Mrot[0][2] = SO->VolPar->ROTATE_MATVEC[2];
1578 Delta[0] = SO->VolPar->ROTATE_MATVEC[3];
1579 Mrot[1][0] = SO->VolPar->ROTATE_MATVEC[4];
1580 Mrot[1][1] = SO->VolPar->ROTATE_MATVEC[5];
1581 Mrot[1][2] = SO->VolPar->ROTATE_MATVEC[6];
1582 Delta[1] = SO->VolPar->ROTATE_MATVEC[7];
1583 Mrot[2][0] = SO->VolPar->ROTATE_MATVEC[8];
1584 Mrot[2][1] = SO->VolPar->ROTATE_MATVEC[9];
1585 Mrot[2][2] = SO->VolPar->ROTATE_MATVEC[10];
1586 Delta[2] = SO->VolPar->ROTATE_MATVEC[11];
1587
1588 NetShift[0] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[0];
1589 NetShift[1] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[1];
1590 NetShift[2] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[2];
1591
1592 for (i=0; i < SO->N_Node; ++i) {
1593 id = ND * i;
1594 /* zero the center */
1595 x = SO->NodeList[id ] - SO->VolPar->ROTATE_CENTER_OLD[0];
1596 y = SO->NodeList[id+1] - SO->VolPar->ROTATE_CENTER_OLD[1];
1597 z = SO->NodeList[id+2] - SO->VolPar->ROTATE_CENTER_OLD[1];
1598
1599 /* Apply the rotation matrix XYZn = Mrot x XYZ*/
1600 SO->NodeList[id ] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
1601 SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
1602 SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
1603
1604 /*apply netshift*/
1605 SO->NodeList[id ] += NetShift[0];
1606 SO->NodeList[id+1] += NetShift[1];
1607 SO->NodeList[id+2] += NetShift[2];
1608 }
1609 SO->ROTATE_APPLIED = YUP;
1610 } else
1611 SO->ROTATE_APPLIED = NOPE;
1612 #endif
1613
1614 SUMA_RETURN (!Bad);
1615 }
1616
1617 /*! The following functions are adaptations from thd_coords.c
1618 They are modified to avoid the use of dset data structure
1619 which is not fully preserved in the Surface Object Data structure.
1620
1621 Stolen Comment:
1622 ====================================================================
1623 3D coordinate conversion routines;
1624 tags for coordinate systems:
1625 fdind = FD_brick voxel indices (ints)
1626 fdfind = FD_brick voxel indices (floats - added 30 Aug 2001)
1627 3dind = THD_3dim_dataset voxel indices (ints)
1628 3dfind = THD_3dim_dataset floating point voxel indices
1629 (used for subvoxel resolution)
1630 3dmm = THD_3dim_dataset millimetric coordinates (floats)
1631 dicomm = DICOM 3.0 millimetric coordinates (floats)
1632
1633 The 3dmm coordinate measurements are oriented the same way
1634 as the dicomm coordinates (which means that the ??del values
1635 can be negative if the voxel axes are reflected), but the 3dmm
1636 coordinates are in general a permutation of the DICOM coordinates.
1637
1638 These routines all take as input and produce as output
1639 THD_?vec3 (i=int,f=float) structs.
1640 ======================================================================
1641
1642 */
1643
SUMA_THD_3dfind_to_3dmm(SUMA_SurfaceObject * SO,THD_fvec3 iv)1644 THD_fvec3 SUMA_THD_3dfind_to_3dmm( SUMA_SurfaceObject *SO,
1645 THD_fvec3 iv )
1646 {
1647 static char FuncName[]={"SUMA_THD_3dfind_to_3dmm"};
1648 THD_fvec3 fv ;
1649
1650 SUMA_ENTRY;
1651
1652 fv.xyz[0] = SO->VolPar->xorg + iv.xyz[0] * SO->VolPar->dx ;
1653 fv.xyz[1] = SO->VolPar->yorg + iv.xyz[1] * SO->VolPar->dy ;
1654 fv.xyz[2] = SO->VolPar->zorg + iv.xyz[2] * SO->VolPar->dz ;
1655 SUMA_RETURN(fv) ;
1656 }
1657
1658 /*!
1659 Same as SUMA_THD_3dfind_to_3dmm, but without needing SO
1660 */
SUMA_THD_3dfind_to_3dmm_vp(SUMA_VOLPAR * vp,THD_fvec3 iv)1661 THD_fvec3 SUMA_THD_3dfind_to_3dmm_vp( SUMA_VOLPAR *vp,
1662 THD_fvec3 iv )
1663 {
1664 static char FuncName[]={"SUMA_THD_3dfind_to_3dmm_vp"};
1665 THD_fvec3 fv ;
1666
1667 SUMA_ENTRY;
1668
1669 fv.xyz[0] = vp->xorg + iv.xyz[0] * vp->dx ;
1670 fv.xyz[1] = vp->yorg + iv.xyz[1] * vp->dy ;
1671 fv.xyz[2] = vp->zorg + iv.xyz[2] * vp->dz ;
1672 SUMA_RETURN(fv) ;
1673 }
1674
1675 /*!------------------------------------------------------------------*/
1676
SUMA_THD_3dind_to_3dmm(SUMA_SurfaceObject * SO,THD_ivec3 iv)1677 THD_fvec3 SUMA_THD_3dind_to_3dmm( SUMA_SurfaceObject *SO,
1678 THD_ivec3 iv )
1679 {
1680 static char FuncName[]={"SUMA_THD_3dind_to_3dmm"};
1681 THD_fvec3 fv ;
1682
1683 SUMA_ENTRY;
1684
1685 fv.xyz[0] = SO->VolPar->xorg + iv.ijk[0] * SO->VolPar->dx ;
1686 fv.xyz[1] = SO->VolPar->yorg + iv.ijk[1] * SO->VolPar->dy ;
1687 fv.xyz[2] = SO->VolPar->zorg + iv.ijk[2] * SO->VolPar->dz ;
1688 SUMA_RETURN(fv) ;
1689 }
1690
1691 /*!--------------------------------------------------------------------*/
1692
SUMA_THD_3dmm_to_3dfind(SUMA_SurfaceObject * SO,THD_fvec3 fv)1693 THD_fvec3 SUMA_THD_3dmm_to_3dfind( SUMA_SurfaceObject *SO ,
1694 THD_fvec3 fv )
1695 {
1696 static char FuncName[]={"SUMA_THD_3dmm_to_3dfind"};
1697 THD_fvec3 iv ;
1698
1699 SUMA_ENTRY;
1700
1701
1702 iv.xyz[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx ;
1703 iv.xyz[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy ;
1704 iv.xyz[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz ;
1705
1706 if( iv.xyz[0] < 0 ) iv.xyz[0] = 0 ;
1707 else if( iv.xyz[0] > SO->VolPar->nx-1 ) iv.xyz[0] = SO->VolPar->nx-1 ;
1708
1709 if( iv.xyz[1] < 0 ) iv.xyz[1] = 0 ;
1710 else if( iv.xyz[1] > SO->VolPar->ny-1 ) iv.xyz[1] = SO->VolPar->ny-1 ;
1711
1712 if( iv.xyz[2] < 0 ) iv.xyz[2] = 0 ;
1713 else if( iv.xyz[2] > SO->VolPar->nz-1 ) iv.xyz[2] = SO->VolPar->nz-1 ;
1714
1715 SUMA_RETURN(iv) ;
1716 }
1717
1718 /*!--------------------------------------------------------------------*/
1719
SUMA_THD_3dmm_to_3dind(SUMA_SurfaceObject * SO,THD_fvec3 fv)1720 THD_ivec3 SUMA_THD_3dmm_to_3dind( SUMA_SurfaceObject *SO ,
1721 THD_fvec3 fv )
1722 {
1723 static char FuncName[]={"SUMA_THD_3dmm_to_3dind"};
1724 THD_ivec3 iv ;
1725
1726 SUMA_ENTRY;
1727
1728 iv.ijk[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx + 0.499 ;
1729 iv.ijk[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy + 0.499 ;
1730 iv.ijk[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz + 0.499 ;
1731
1732 if( iv.ijk[0] < 0 ) iv.ijk[0] = 0 ;
1733 else if( iv.ijk[0] > SO->VolPar->nx-1 ) iv.ijk[0] = SO->VolPar->nx-1 ;
1734
1735 if( iv.ijk[1] < 0 ) iv.ijk[1] = 0 ;
1736 else if( iv.ijk[1] > SO->VolPar->ny-1 ) iv.ijk[1] = SO->VolPar->ny-1 ;
1737
1738 if( iv.ijk[2] < 0 ) iv.ijk[2] = 0 ;
1739 else if( iv.ijk[2] > SO->VolPar->nz-1 ) iv.ijk[2] = SO->VolPar->nz-1 ;
1740
1741 SUMA_RETURN(iv) ;
1742 }
1743
SUMA_THD_3dmm_to_3dind_warn(SUMA_SurfaceObject * SO,THD_fvec3 fv,int * out)1744 THD_ivec3 SUMA_THD_3dmm_to_3dind_warn( SUMA_SurfaceObject *SO ,
1745 THD_fvec3 fv , int *out)
1746 {
1747 static char FuncName[]={"SUMA_THD_3dmm_to_3dind_warn"};
1748 THD_ivec3 iv ;
1749
1750 SUMA_ENTRY;
1751 *out = 0;
1752
1753 iv.ijk[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx + 0.499 ;
1754 iv.ijk[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy + 0.499 ;
1755 iv.ijk[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz + 0.499 ;
1756
1757 if( iv.ijk[0] < 0 ) { iv.ijk[0] = 0 ; *out = 1; }
1758 else if( iv.ijk[0] > SO->VolPar->nx-1 ) { iv.ijk[0] = SO->VolPar->nx-1 ; *out = 1; }
1759
1760 if( iv.ijk[1] < 0 ) { iv.ijk[1] = 0 ; *out = 1; }
1761 else if( iv.ijk[1] > SO->VolPar->ny-1 ) { iv.ijk[1] = SO->VolPar->ny-1 ; *out = 1; }
1762
1763 if( iv.ijk[2] < 0 ) { iv.ijk[2] = 0 ; *out = 1; }
1764 else if( iv.ijk[2] > SO->VolPar->nz-1 ) { iv.ijk[2] = SO->VolPar->nz-1 ; *out = 1; }
1765
1766 SUMA_RETURN(iv) ;
1767 }
1768 /*!
1769 \brief how many voxels in each of the RL AP IS directions
1770 */
SUMA_VolDims(THD_3dim_dataset * dset,int * nRL,int * nAP,int * nIS)1771 void SUMA_VolDims(THD_3dim_dataset *dset, int *nRL, int *nAP, int *nIS)
1772 {
1773 static char FuncName[]={"SUMA_VolDims"};
1774
1775 SUMA_ENTRY;
1776
1777 *nRL = *nAP = *nIS = -1;
1778
1779 if (!dset) {
1780 SUMA_SL_Err("NULL dset");
1781 SUMA_RETURNe;
1782 }
1783
1784 switch( dset->daxes->xxorient ){
1785 case ORI_R2L_TYPE:
1786 case ORI_L2R_TYPE: *nRL = DSET_NX(dset) ; break ;
1787 case ORI_P2A_TYPE:
1788 case ORI_A2P_TYPE: *nAP = DSET_NX(dset) ; break ;
1789 case ORI_I2S_TYPE:
1790 case ORI_S2I_TYPE: *nIS = DSET_NX(dset) ; break ;
1791 }
1792
1793 switch( dset->daxes->yyorient ){
1794 case ORI_R2L_TYPE:
1795 case ORI_L2R_TYPE: *nRL = DSET_NY(dset) ; break ;
1796 case ORI_P2A_TYPE:
1797 case ORI_A2P_TYPE: *nAP = DSET_NY(dset) ; break ;
1798 case ORI_I2S_TYPE:
1799 case ORI_S2I_TYPE: *nIS = DSET_NY(dset) ; break ;
1800 }
1801
1802 switch( dset->daxes->zzorient ){
1803 case ORI_R2L_TYPE:
1804 case ORI_L2R_TYPE: *nRL = DSET_NZ(dset) ; break ;
1805 case ORI_P2A_TYPE:
1806 case ORI_A2P_TYPE: *nAP = DSET_NZ(dset) ; break ;
1807 case ORI_I2S_TYPE:
1808 case ORI_S2I_TYPE: *nIS = DSET_NZ(dset) ; break ;
1809 }
1810
1811 SUMA_RETURNe;
1812 }
1813
1814 /*!---------------------------------------------------------------------
1815 convert from input image oriented x,y,z to Dicom x,y,z
1816 (x axis = R->L , y axis = A->P , z axis = I->S)
1817
1818 N.B.: image distances are oriented the same as Dicom,
1819 just in a permuted order.
1820 -----------------------------------------------------------------------*/
SUMA_THD_3dmm_to_dicomm(int xxorient,int yyorient,int zzorient,THD_fvec3 imv)1821 THD_fvec3 SUMA_THD_3dmm_to_dicomm( int xxorient, int yyorient, int zzorient,
1822 THD_fvec3 imv )
1823 {
1824 static char FuncName[]={"SUMA_THD_3dmm_to_dicomm"};
1825 THD_fvec3 dicv ;
1826 float xim,yim,zim , xdic = 0.0, ydic = 0.0, zdic = 0.0 ;
1827
1828 SUMA_ENTRY;
1829
1830 xim = imv.xyz[0] ; yim = imv.xyz[1] ; zim = imv.xyz[2] ;
1831
1832 switch( xxorient ){
1833 case ORI_R2L_TYPE:
1834 case ORI_L2R_TYPE: xdic = xim ; break ;
1835 case ORI_P2A_TYPE:
1836 case ORI_A2P_TYPE: ydic = xim ; break ;
1837 case ORI_I2S_TYPE:
1838 case ORI_S2I_TYPE: zdic = xim ; break ;
1839
1840 default:
1841 fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
1842 exit (1);
1843 }
1844
1845 switch( yyorient ){
1846 case ORI_R2L_TYPE:
1847 case ORI_L2R_TYPE: xdic = yim ; break ;
1848 case ORI_P2A_TYPE:
1849 case ORI_A2P_TYPE: ydic = yim ; break ;
1850 case ORI_I2S_TYPE:
1851 case ORI_S2I_TYPE: zdic = yim ; break ;
1852
1853 default:
1854 fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
1855 exit (1);
1856 }
1857
1858 switch( zzorient ){
1859 case ORI_R2L_TYPE:
1860 case ORI_L2R_TYPE: xdic = zim ; break ;
1861 case ORI_P2A_TYPE:
1862 case ORI_A2P_TYPE: ydic = zim ; break ;
1863 case ORI_I2S_TYPE:
1864 case ORI_S2I_TYPE: zdic = zim ; break ;
1865
1866 default:
1867 fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
1868 exit (1);
1869 }
1870
1871 dicv.xyz[0] = xdic ; dicv.xyz[1] = ydic ; dicv.xyz[2] = zdic ;
1872 SUMA_RETURN(dicv) ;
1873 }
1874
1875 /*!---------------------------------------------------------------------
1876 convert to input image oriented x,y,z from Dicom x,y,z
1877 -----------------------------------------------------------------------*/
1878
SUMA_THD_dicomm_to_3dmm(int xxorient,int yyorient,int zzorient,THD_fvec3 dicv)1879 THD_fvec3 SUMA_THD_dicomm_to_3dmm( int xxorient, int yyorient, int zzorient ,
1880 THD_fvec3 dicv )
1881 {
1882 static char FuncName[]={"SUMA_THD_dicomm_to_3dmm"};
1883 THD_fvec3 imv ;
1884 float xim,yim,zim , xdic,ydic,zdic ;
1885
1886 SUMA_ENTRY;
1887
1888 xdic = dicv.xyz[0] ; ydic = dicv.xyz[1] ; zdic = dicv.xyz[2] ;
1889
1890 switch( xxorient ){
1891 case ORI_R2L_TYPE:
1892 case ORI_L2R_TYPE: xim = xdic ; break ;
1893 case ORI_P2A_TYPE:
1894 case ORI_A2P_TYPE: xim = ydic ; break ;
1895 case ORI_I2S_TYPE:
1896 case ORI_S2I_TYPE: xim = zdic ; break ;
1897
1898 default:
1899 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
1900 exit (1);
1901 }
1902
1903 switch( yyorient ){
1904 case ORI_R2L_TYPE:
1905 case ORI_L2R_TYPE: yim = xdic ; break ;
1906 case ORI_P2A_TYPE:
1907 case ORI_A2P_TYPE: yim = ydic ; break ;
1908 case ORI_I2S_TYPE:
1909 case ORI_S2I_TYPE: yim = zdic ; break ;
1910
1911 default:
1912 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
1913 exit (1);
1914 }
1915
1916 switch( zzorient ){
1917 case ORI_R2L_TYPE:
1918 case ORI_L2R_TYPE: zim = xdic ; break ;
1919 case ORI_P2A_TYPE:
1920 case ORI_A2P_TYPE: zim = ydic ; break ;
1921 case ORI_I2S_TYPE:
1922 case ORI_S2I_TYPE: zim = zdic ; break ;
1923
1924 default:
1925 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
1926 exit (1);
1927 }
1928
1929 imv.xyz[0] = xim ; imv.xyz[1] = yim ; imv.xyz[2] = zim ;
1930 SUMA_RETURN(imv) ;
1931 }
1932 /*!
1933 \brief changes orientation codes to string
1934 classic RAI codes: ORI_R2L_TYPE, ORI_A2P_TYPE, ORI_I2S_TYPE
1935 would result in orstr being set to: RAILPS
1936 */
SUMA_orcode_to_orstring(int xxorient,int yyorient,int zzorient,char * orstr)1937 void SUMA_orcode_to_orstring (int xxorient,int yyorient,int zzorient, char *orstr)
1938 {
1939 static char FuncName[]={"SUMA_orcode_to_orstring"};
1940
1941 SUMA_ENTRY;
1942
1943 if (!orstr) { SUMA_SL_Err("NULL string"); SUMA_RETURNe; }
1944
1945 orstr[0]='\0';
1946 switch( xxorient ){
1947 case ORI_R2L_TYPE: orstr[0] = 'R'; orstr[3] = 'L'; break ;
1948 case ORI_L2R_TYPE: orstr[0] = 'L'; orstr[3] = 'R'; break ;
1949 case ORI_P2A_TYPE: orstr[0] = 'P'; orstr[3] = 'A'; break ;
1950 case ORI_A2P_TYPE: orstr[0] = 'A'; orstr[3] = 'P'; break ;
1951 case ORI_I2S_TYPE: orstr[0] = 'I'; orstr[3] = 'S'; break ;
1952 case ORI_S2I_TYPE: orstr[0] = 'S'; orstr[3] = 'I'; break ;
1953
1954 default:
1955 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n") ;
1956 SUMA_RETURNe;
1957 }
1958
1959 switch( yyorient ){
1960 case ORI_R2L_TYPE: orstr[1] = 'R'; orstr[4] = 'L'; break ;
1961 case ORI_L2R_TYPE: orstr[1] = 'L'; orstr[4] = 'R'; break ;
1962 case ORI_P2A_TYPE: orstr[1] = 'P'; orstr[4] = 'A'; break ;
1963 case ORI_A2P_TYPE: orstr[1] = 'A'; orstr[4] = 'P'; break ;
1964 case ORI_I2S_TYPE: orstr[1] = 'I'; orstr[4] = 'S'; break ;
1965 case ORI_S2I_TYPE: orstr[1] = 'S'; orstr[4] = 'I'; break ;
1966
1967 default:
1968 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal yyorient code.\n ") ;
1969 SUMA_RETURNe;
1970 }
1971
1972 switch( zzorient ){
1973 case ORI_R2L_TYPE: orstr[2] = 'R'; orstr[5] = 'L'; break ;
1974 case ORI_L2R_TYPE: orstr[2] = 'L'; orstr[5] = 'R'; break ;
1975 case ORI_P2A_TYPE: orstr[2] = 'P'; orstr[5] = 'A'; break ;
1976 case ORI_A2P_TYPE: orstr[2] = 'A'; orstr[5] = 'P'; break ;
1977 case ORI_I2S_TYPE: orstr[2] = 'I'; orstr[5] = 'S'; break ;
1978 case ORI_S2I_TYPE: orstr[2] = 'S'; orstr[5] = 'I'; break ;
1979
1980 default:
1981 fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal zzorient code.\n ") ;
1982 SUMA_RETURNe;
1983 }
1984
1985 SUMA_RETURNe;
1986 }
SUMA_orstring_to_orcode(char * orstr,int * orient)1987 SUMA_Boolean SUMA_orstring_to_orcode (char *orstr, int *orient)
1988 {
1989 static char FuncName[]={"SUMA_orstring_to_orcode"};
1990 int i;
1991
1992 SUMA_ENTRY;
1993
1994 if (!orstr) { SUMA_SL_Err("NULL string"); SUMA_RETURN(NOPE); }
1995 if (!SUMA_ok_orstring(orstr)) { SUMA_SL_Err("Bad orientation string"); SUMA_RETURN(NOPE); }
1996 for (i=0; i<3; ++i) {
1997 switch (orstr[i]) {
1998 case 'R': orient[i] = ORI_R2L_TYPE; break;
1999 case 'L': orient[i] = ORI_L2R_TYPE; break;
2000 case 'A': orient[i] = ORI_A2P_TYPE; break;
2001 case 'P': orient[i] = ORI_P2A_TYPE; break;
2002 case 'I': orient[i] = ORI_I2S_TYPE; break;
2003 case 'S': orient[i] = ORI_S2I_TYPE; break;
2004 default: fprintf(SUMA_STDERR, " SUMA_orstring_to_orcode: Bad to the bones\n"); SUMA_RETURN(NOPE);
2005 }
2006 }
2007
2008 SUMA_RETURN(YUP);
2009 }
2010
SUMA_ok_orstring(char * orstr)2011 int SUMA_ok_orstring (char *orstr)
2012 {
2013 static char FuncName[]={"SUMA_ok_orstring"};
2014 int i, d[3];
2015
2016 SUMA_ENTRY;
2017
2018 if (!orstr) { SUMA_RETURN(NOPE); }
2019 d[0] = d[1] = d[2] = 0;
2020 for (i=0; i<3; ++i) {
2021 switch (orstr[i]) {
2022 case 'R': ++(d[0]); break;
2023 case 'L': ++(d[0]); break;
2024 case 'A': ++(d[1]); break;
2025 case 'P': ++(d[1]); break;
2026 case 'I': ++(d[2]); break;
2027 case 'S': ++(d[2]); break;
2028 default: SUMA_RETURN(NOPE);
2029 }
2030 }
2031 if (d[0] != 1 || d[1] != 1 || d[2] != 1) SUMA_RETURN(NOPE);
2032
2033 SUMA_RETURN(YUP);
2034 }
SUMA_flip_orient(int xxorient)2035 int SUMA_flip_orient(int xxorient)
2036 {
2037 static char FuncName[]={"SUMA_flip_orient"};
2038
2039 SUMA_ENTRY;
2040
2041 switch( xxorient ){
2042 case ORI_R2L_TYPE: SUMA_RETURN(ORI_L2R_TYPE); break ;
2043 case ORI_L2R_TYPE: SUMA_RETURN(ORI_R2L_TYPE); break ;
2044
2045 case ORI_P2A_TYPE: SUMA_RETURN(ORI_A2P_TYPE); break ;
2046 case ORI_A2P_TYPE: SUMA_RETURN(ORI_P2A_TYPE); break ;
2047
2048 case ORI_I2S_TYPE: SUMA_RETURN(ORI_S2I_TYPE); break ;
2049 case ORI_S2I_TYPE: SUMA_RETURN(ORI_I2S_TYPE); break ;
2050
2051 default:
2052 fprintf(SUMA_STDERR, "SUMA_opposite_orient: illegal zzorient code.\n ") ;
2053 SUMA_RETURN(-1);
2054 }
2055
2056 SUMA_RETURN(-1);
2057
2058 }
2059
SUMA_CoordChange(char * orc_in,char * orc_out,float * XYZ,int N_xyz)2060 SUMA_Boolean SUMA_CoordChange (char *orc_in, char *orc_out, float *XYZ, int N_xyz)
2061 {
2062 static char FuncName[]={"SUMA_CoordChange"};
2063 int i, or_in[3], or_out[3], j, map[3], sgn[3], i3;
2064 float xyz[3];
2065
2066 SUMA_ENTRY;
2067
2068 if (!SUMA_orstring_to_orcode(orc_in, or_in)) {
2069 SUMA_SL_Err("Bad in code");
2070 SUMA_RETURN(NOPE);
2071 }
2072 if (!SUMA_orstring_to_orcode(orc_out, or_out)) {
2073 SUMA_SL_Err("Bad out code");
2074 SUMA_RETURN(NOPE);
2075 }
2076
2077 /* figure out the mapping */
2078 for (j=0; j<3; ++j) {
2079 i = 0;
2080 while (i<3) {
2081 if (or_in[i] == or_out[j] || or_in[i] == SUMA_flip_orient(or_out[j])) {
2082 map[j] = i;
2083 if (or_in[i] == SUMA_flip_orient(or_out[j])) sgn[j] = -1;
2084 else sgn[j] = 1;
2085 i = 3; /* break */
2086 }
2087 ++i;
2088 }
2089 }
2090
2091 for (i=0; i<N_xyz; ++i) {
2092 i3 = 3*i;
2093 xyz[0] = XYZ[i3]; xyz[1] = XYZ[i3+1]; xyz[2] = XYZ[i3+2];
2094 XYZ[i3 ] = sgn[0] * xyz[map[0]];
2095 XYZ[i3+1] = sgn[1] * xyz[map[1]];
2096 XYZ[i3+2] = sgn[2] * xyz[map[2]];
2097 }
2098
2099 SUMA_RETURN(YUP);
2100 }
2101 /*!
2102 \brief takes the origin as entered to to3d and
2103 changes them to the origin field for AFNI header
2104 Based on info in README.attributes and function T3D_save_file_CB
2105 int to3d.c
2106 */
SUMA_originto3d_2_originHEAD(THD_ivec3 orient,THD_fvec3 * origin)2107 void SUMA_originto3d_2_originHEAD(THD_ivec3 orient, THD_fvec3 *origin)
2108 {
2109 static char FuncName[]={"SUMA_originto3d_2_originHEAD"};
2110
2111 SUMA_ENTRY;
2112
2113 origin->xyz[0] = (ORIENT_sign[orient.ijk[0]] == '+')
2114 ? (-origin->xyz[0]) : ( origin->xyz[0]) ;
2115
2116 origin->xyz[1] = (ORIENT_sign[orient.ijk[1]] == '+')
2117 ? (-origin->xyz[1]) : ( origin->xyz[1]) ;
2118
2119 origin->xyz[2] = (ORIENT_sign[orient.ijk[2]] == '+')
2120 ? (-origin->xyz[2]) : ( origin->xyz[2]) ;
2121
2122 SUMA_RETURNe;
2123
2124 }
2125 /*!
2126 \brief takes the size as entered to to3d and
2127 changes them to the delta field for AFNI header
2128 Based on info in README.attributes and function T3D_save_file_CB
2129 int to3d.c
2130 */
2131
SUMA_sizeto3d_2_deltaHEAD(THD_ivec3 orient,THD_fvec3 * delta)2132 void SUMA_sizeto3d_2_deltaHEAD(THD_ivec3 orient, THD_fvec3 *delta)
2133 {
2134 static char FuncName[]={"SUMA_sizeto3d_2_deltaHEAD"};
2135
2136 SUMA_ENTRY;
2137
2138 delta->xyz[0] = (ORIENT_sign[orient.ijk[0]] == '+')
2139 ? (delta->xyz[0]) : ( -delta->xyz[0]) ;
2140
2141 delta->xyz[1] = (ORIENT_sign[orient.ijk[1]] == '+')
2142 ? (delta->xyz[1]) : ( -delta->xyz[1]) ;
2143
2144 delta->xyz[2] = (ORIENT_sign[orient.ijk[2]] == '+')
2145 ? (delta->xyz[2]) : ( -delta->xyz[2]) ;
2146
2147 SUMA_RETURNe;
2148
2149 }
2150
2151 /*!
2152 \brief SUMA_Boolean SUMA_vec_3dfind_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
2153 SUMA_Boolean SUMA_vec_3dmm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
2154 SUMA_Boolean SUMA_vec_dicomm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
2155 SUMA_Boolean SUMA_vec_3dfind_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
2156 SUMA_Boolean SUMA_vec_3dmm_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
2157 SUMA_Boolean SUMA_vec_dicomm_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
2158 a set of functions to change coordinate systems.
2159
2160 \param NodeList (float *) vector of consecutive xyz (N_Node x 3) values
2161 \param N_Node (int) number of xyz triplets in NodeList
2162 \param VolPar (SUMA_VOLPAR *) volume parent structure (output of SUMA_VolPar_Attr )
2163
2164 3dfind = float voxel/node index coordinate
2165 3dmm = float voxel/node coordinate in volume space
2166 dicomm = float voxel/node coordinate in RAI space
2167
2168 see also SUMA_THD_ functions ...
2169 */
SUMA_vec_3dfind_to_3dmm(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2170 SUMA_Boolean SUMA_vec_3dfind_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
2171 {
2172 static char FuncName[]={"SUMA_vec_3dfind_to_3dmm"};
2173 THD_fvec3 fv, iv;
2174 int i, id;
2175 SUMA_SurfaceObject SO;
2176
2177 SUMA_ENTRY;
2178
2179 if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2180 /* create dummy struct */
2181 SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
2182
2183 for (i=0; i < SO.N_Node; ++i) {
2184 id = i * SO.NodeDim;
2185 /* change float indices to mm coords */
2186 iv.xyz[0] = SO.NodeList[id] ;
2187 iv.xyz[1] = SO.NodeList[id+1] ;
2188 iv.xyz[2] = SO.NodeList[id+2] ;
2189 fv = SUMA_THD_3dfind_to_3dmm( &SO, iv );
2190 SO.NodeList[id] = fv.xyz[0];
2191 SO.NodeList[id+1] = fv.xyz[1];
2192 SO.NodeList[id+2] = fv.xyz[2];
2193 }
2194
2195 SUMA_RETURN(YUP);
2196 }
2197
SUMA_vec_3dmm_to_3dfind(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2198 SUMA_Boolean SUMA_vec_3dmm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
2199 {
2200 static char FuncName[]={"SUMA_vec_3dmm_to_3dfind"};
2201 THD_fvec3 fv, iv;
2202 int i, id;
2203 SUMA_SurfaceObject SO;
2204
2205 SUMA_ENTRY;
2206
2207 if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2208 /* create dummy struct */
2209 SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
2210
2211 for (i=0; i < SO.N_Node; ++i) {
2212 id = i * SO.NodeDim;
2213 /* change float indices to mm coords */
2214 iv.xyz[0] = SO.NodeList[id] ;
2215 iv.xyz[1] = SO.NodeList[id+1] ;
2216 iv.xyz[2] = SO.NodeList[id+2] ;
2217 fv = SUMA_THD_3dmm_to_3dfind( &SO, iv );
2218 /* fprintf(SUMA_STDERR,"%s: In[%f %f %f] Out[%f %f %f]\n",
2219 FuncName, iv.xyz[0], iv.xyz[1], iv.xyz[2],
2220 fv.xyz[0], fv.xyz[1], fv.xyz[2]); */
2221 SO.NodeList[id] = fv.xyz[0];
2222 SO.NodeList[id+1] = fv.xyz[1];
2223 SO.NodeList[id+2] = fv.xyz[2];
2224 }
2225
2226 SUMA_RETURN(YUP);
2227 }
2228
SUMA_vec_dicomm_to_3dfind(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2229 SUMA_Boolean SUMA_vec_dicomm_to_3dfind (float *NodeList,
2230 int N_Node, SUMA_VOLPAR *VolPar)
2231 {
2232 static char FuncName[]={"SUMA_vec_dicomm_to_3dfind"};
2233
2234 SUMA_ENTRY;
2235
2236 if (!NodeList || !VolPar) {
2237 SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2238
2239
2240 if (!SUMA_vec_dicomm_to_3dmm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
2241 if (!SUMA_vec_3dmm_to_3dfind(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
2242
2243 SUMA_RETURN(YUP);
2244 }
2245
SUMA_vec_3dfind_to_dicomm(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2246 SUMA_Boolean SUMA_vec_3dfind_to_dicomm (float *NodeList, int N_Node,
2247 SUMA_VOLPAR *VolPar)
2248 {
2249 static char FuncName[]={"SUMA_vec_3dfind_to_dicomm"};
2250
2251 SUMA_ENTRY;
2252
2253 if (!NodeList || !VolPar) {
2254 SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2255
2256 if (!SUMA_vec_3dfind_to_3dmm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
2257 if (!SUMA_vec_3dmm_to_dicomm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
2258
2259 SUMA_RETURN(YUP);
2260 }
2261
SUMA_vec_3dmm_to_dicomm(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2262 SUMA_Boolean SUMA_vec_3dmm_to_dicomm (float *NodeList, int N_Node,
2263 SUMA_VOLPAR *VolPar)
2264 {
2265 static char FuncName[]={"SUMA_vec_3dmm_to_dicomm"};
2266 THD_fvec3 fv, iv;
2267 int i, id, NodeDim=3;
2268
2269 SUMA_ENTRY;
2270
2271 if (!NodeList || !VolPar) {
2272 SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2273
2274 for (i=0; i < N_Node; ++i) {
2275 id = i * NodeDim;
2276 iv.xyz[0] = NodeList[id] ;
2277 iv.xyz[1] = NodeList[id+1] ;
2278 iv.xyz[2] = NodeList[id+2] ;
2279
2280 /* change mm to RAI coords */
2281 fv = SUMA_THD_3dmm_to_dicomm( VolPar->xxorient,
2282 VolPar->yyorient,
2283 VolPar->zzorient, iv );
2284 NodeList[id] = fv.xyz[0];
2285 NodeList[id+1] = fv.xyz[1];
2286 NodeList[id+2] = fv.xyz[2];
2287 }
2288
2289 SUMA_RETURN(YUP);
2290 }
2291
SUMA_vec_dicomm_to_3dmm(float * NodeList,int N_Node,SUMA_VOLPAR * VolPar)2292 SUMA_Boolean SUMA_vec_dicomm_to_3dmm (float *NodeList, int N_Node,
2293 SUMA_VOLPAR *VolPar)
2294 {
2295 static char FuncName[]={"SUMA_vec_dicomm_to_3dmm"};
2296 THD_fvec3 fv, iv;
2297 int i, id, NodeDim=3;
2298
2299 SUMA_ENTRY;
2300
2301 if (!NodeList || !VolPar) {
2302 SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
2303
2304 for (i=0; i < N_Node; ++i) {
2305 id = i * NodeDim;
2306
2307 iv.xyz[0] = NodeList[id] ;
2308 iv.xyz[1] = NodeList[id+1] ;
2309 iv.xyz[2] = NodeList[id+2] ;
2310
2311 /* change mm to RAI coords */
2312 fv = SUMA_THD_dicomm_to_3dmm( VolPar->xxorient,
2313 VolPar->yyorient, VolPar->zzorient,iv );
2314 /* fprintf(SUMA_STDERR,"%s: In[%f %f %f] Out[%f %f %f]\n",
2315 FuncName, iv.xyz[0], iv.xyz[1], iv.xyz[2],
2316 fv.xyz[0], fv.xyz[1], fv.xyz[2]); */
2317 NodeList[id] = fv.xyz[0];
2318 NodeList[id+1] = fv.xyz[1];
2319 NodeList[id+2] = fv.xyz[2];
2320 }
2321
2322 SUMA_RETURN(YUP);
2323 }
2324
SUMA_THD_3dfind_to_dicomm(THD_3dim_dataset * dset,float ii,float jj,float kk,float * xyz)2325 SUMA_Boolean SUMA_THD_3dfind_to_dicomm(THD_3dim_dataset *dset,
2326 float ii, float jj, float kk,
2327 float *xyz)
2328 {
2329 THD_fvec3 fv0, fv;
2330 fv0.xyz[0]=ii;
2331 fv0.xyz[1]=jj;
2332 fv0.xyz[2]=kk;
2333 fv = THD_3dfind_to_3dmm(dset, fv0);
2334 fv0 = THD_3dmm_to_dicomm(dset, fv);
2335 xyz[0] = fv0.xyz[0];
2336 xyz[1] = fv0.xyz[1];
2337 xyz[2] = fv0.xyz[2];
2338 return(1);
2339 }
2340
SUMA_THD_dicomm_to_3dfind(THD_3dim_dataset * dset,float RR,float AA,float II,float * ijk)2341 SUMA_Boolean SUMA_THD_dicomm_to_3dfind(THD_3dim_dataset *dset,
2342 float RR, float AA, float II,
2343 float *ijk)
2344 {
2345 THD_fvec3 fv0, fv;
2346 fv0.xyz[0]=RR;
2347 fv0.xyz[1]=AA;
2348 fv0.xyz[2]=II;
2349 fv = THD_dicomm_to_3dmm(dset, fv0);
2350 fv0 = THD_3dmm_to_3dfind(dset, fv);
2351 ijk[0] = fv0.xyz[0];
2352 ijk[1] = fv0.xyz[1];
2353 ijk[2] = fv0.xyz[2];
2354 return(1);
2355 }
2356
SUMA_THD_dicomm_to_1dind(THD_3dim_dataset * dset,float RR,float AA,float II,int * ijk)2357 int SUMA_THD_dicomm_to_1dind(THD_3dim_dataset *dset,
2358 float RR, float AA, float II,
2359 int *ijk)
2360 {
2361 THD_fvec3 fv0, fv;
2362 fv0.xyz[0]=RR;
2363 fv0.xyz[1]=AA;
2364 fv0.xyz[2]=II;
2365 fv = THD_dicomm_to_3dmm(dset, fv0);
2366 fv0 = THD_3dmm_to_3dfind(dset, fv);
2367 if (ijk) {
2368 ijk[0] = fv0.xyz[0];
2369 ijk[1] = fv0.xyz[1];
2370 ijk[2] = fv0.xyz[2];
2371 }
2372 return((int)fv0.xyz[0]+(int)fv0.xyz[1]*DSET_NX(dset)+
2373 (int)fv0.xyz[2]*DSET_NX(dset)*DSET_NY(dset));
2374 }
2375
2376 /*!
2377 \brief transforms XYZ coordinates from AFNI'S RAI
2378 talairach space to MNI space in in RAI or (LPI). see, e.g.:
2379 https://imaging.mrc-cbu.cam.ac.uk/imaging/MniTalairach
2380 Also see Bob's THD_tta_to_mni in thd_mnicoords.c
2381
2382 \param NodeList (float *) vector of coordinates, 3 * N_Node long
2383 \param N_Node (int) number of nodes in vector above
2384 \param Coord (char *) one of two options
2385 "RAI" meaning the MNI coordinates are to be in RAI
2386 "LPI" meaning the MNI coordinates are to be in LPI
2387 \return flag (SUMA_Boolean) NOPE = failed
2388
2389 - This function should be rewritten to call THD_tta_to_mni (although you'll have to deal
2390 with the flipping option
2391
2392
2393 */
SUMA_AFNItlrc_toMNI(float * NodeList,int N_Node,char * Coord)2394 SUMA_Boolean SUMA_AFNItlrc_toMNI(float *NodeList, int N_Node, char *Coord)
2395 {
2396 static char FuncName[]={"SUMA_AFNItlrc_toMNI"};
2397 SUMA_Boolean DoFlip = NOPE;
2398 int i, i3;
2399 float mx = 0.0,my = 0.0,mz = 0.0, tx = 0.0,ty = 0.0,tz = 0.0 ;
2400 SUMA_Boolean LocalHead = NOPE;
2401
2402 SUMA_ENTRY;
2403
2404 if (strcmp(Coord,"RAI") == 0) DoFlip = NOPE;
2405 else if (strcmp(Coord,"LPI") == 0) DoFlip = YUP;
2406 else {
2407 SUMA_S_Err("Can't do. Either RAI or LPI");
2408 SUMA_RETURN(NOPE);
2409 }
2410
2411
2412 for (i=0; i< N_Node; ++i) {
2413 i3 = 3*i;
2414 if (DoFlip) {
2415 if (!i) SUMA_LH("Flipping to LPI");
2416 tx = - NodeList[i3]; ty = -NodeList[i3+1] ; /* flip xy from RAI to LPI */
2417 tz = NodeList[i3+2] ;
2418 } else {
2419 if (!i) SUMA_LH("No Flipping, RAI maintained");
2420 tx = NodeList[i3]; ty = NodeList[i3+1] ; /* flip xy from RAI to LPI */
2421 tz = NodeList[i3+2] ;
2422 }
2423 mx = 1.01010 * tx ;
2424 my = 1.02962 * ty - 0.05154 * tz ;
2425 mz = 0.05434 * ty + 1.08554 * tz ;
2426 if( mz < 0.0 ) mz *= 1.09523 ;
2427 NodeList[i3] = mx;
2428 NodeList[i3+1] = my;
2429 NodeList[i3+2] = mz;
2430 }
2431
2432
2433 SUMA_RETURN(YUP);
2434 }
2435 /*!
2436
2437 \brief transforms XYZ coordinates by transfrom in warp.
2438 ans = SUMA_AFNI_forward_warp_xyz( warp , XYZv, N);
2439
2440 \param warp (THD_warp *) afni warp structure
2441 \param XYZv (float *) pointer to coordinates vector.
2442 \param N (int) number of XYZ triplets in XYZv.
2443 The ith X,Y, Z coordinates would be XYZv[3i], XYZv[3i+1],XYZv[3i+2]
2444 \return ans (YUP/NOPE) NOPE: NULL warp or bad warp->type
2445
2446 - The following functions are adapted from afni.c & Vecwarp.c
2447 */
SUMA_AFNI_forward_warp_xyz(THD_warp * warp,float * xyzv,int N)2448 SUMA_Boolean SUMA_AFNI_forward_warp_xyz( THD_warp * warp , float *xyzv, int N)
2449 {
2450 static char FuncName[]={"SUMA_AFNI_forward_warp_xyz"};
2451 THD_fvec3 new_fv, old_fv;
2452 int i, i3;
2453 SUMA_Boolean LocalHead = NOPE;
2454
2455 SUMA_ENTRY;
2456
2457 if( warp == NULL ) {
2458 fprintf (SUMA_STDERR, "Error %s: NULL warp.\n", FuncName);
2459 SUMA_RETURN (NOPE) ;
2460 }
2461
2462 switch( warp->type ){
2463
2464 default:
2465 fprintf (SUMA_STDERR, "Error %s: bad warp->type\n", FuncName);
2466 SUMA_RETURN (NOPE);
2467 break ;
2468
2469 case WARP_TALAIRACH_12_TYPE:{
2470 THD_linear_mapping map ;
2471 int iw ;
2472
2473 /* forward transform each possible case,
2474 and test if result is in bot..top of defined map */
2475
2476 for (i=0; i < N; ++i) {
2477 i3 = 3*i;
2478 old_fv.xyz[0] = xyzv[i3];
2479 old_fv.xyz[1] = xyzv[i3+1];
2480 old_fv.xyz[2] = xyzv[i3+2];
2481
2482 for( iw=0 ; iw < 12 ; iw++ ){
2483 map = warp->tal_12.warp[iw] ;
2484 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
2485
2486 if( new_fv.xyz[0] >= map.bot.xyz[0] &&
2487 new_fv.xyz[1] >= map.bot.xyz[1] &&
2488 new_fv.xyz[2] >= map.bot.xyz[2] &&
2489 new_fv.xyz[0] <= map.top.xyz[0] &&
2490 new_fv.xyz[1] <= map.top.xyz[1] &&
2491 new_fv.xyz[2] <= map.top.xyz[2] ) break ; /* leave loop */
2492 }
2493 xyzv[i3] = new_fv.xyz[0];
2494 xyzv[i3+1] = new_fv.xyz[1];
2495 xyzv[i3+2] = new_fv.xyz[2];
2496
2497 }
2498 }
2499 break ;
2500
2501 case WARP_AFFINE_TYPE:{
2502 THD_linear_mapping map = warp->rig_bod.warp ;
2503 SUMA_LH("Have WarpAffineType");
2504 if (LocalHead) {
2505 for (i=0; i < 3; ++i) fprintf(SUMA_STDERR,"%.5f %.5f %.5f %.5f\n",
2506 map.mfor.mat[i][0], map.mfor.mat[i][1], map.mfor.mat[i][2], map.bvec.xyz[i]);
2507 }
2508 for (i=0; i < N; ++i) {
2509 i3 = 3*i;
2510 old_fv.xyz[0] = xyzv[i3];
2511 old_fv.xyz[1] = xyzv[i3+1];
2512 old_fv.xyz[2] = xyzv[i3+2];
2513 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
2514 xyzv[i3] = new_fv.xyz[0];
2515 xyzv[i3+1] = new_fv.xyz[1];
2516 xyzv[i3+2] = new_fv.xyz[2];
2517 }
2518 }
2519 break ;
2520
2521 }
2522 SUMA_RETURN(YUP);
2523 }
2524