1 #include "mrilib.h"
2
3 /*--------------------------------------------------------------------------
4 Program to transform 3-vectors based on warps stored in AFNI .HEAD files.
5 Transformations are stored in +view.HEAD files (view=acpc or tlrc) between
6 the +orig and +view coordinate systems.
7
8 24 Oct 2001: creation by RWCox
9 ----------------------------------------------------------------------------*/
10
11 /*** Some prototypes ***/
12
13 THD_fvec3 myAFNI_forward_warp_vector ( THD_warp * , THD_fvec3 ) ;
14 THD_fvec3 myAFNI_backward_warp_vector( THD_warp * , THD_fvec3 ) ;
15
16 #if 0
17 static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *, THD_fvec3 ) ;
18 static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *, THD_fvec3 ) ;
19 #endif
20
21 #undef DEBUG
22
23 /*--------------------------------------------------------------------------*/
24
Syntax(void)25 static void Syntax(void)
26 {
27 printf(
28 "Usage: Vecwarp [options]\n"
29 "Transforms (warps) a list of 3-vectors into another list of 3-vectors\n"
30 "according to the options. Error messages, warnings, and informational\n"
31 "messages are written to stderr. If a fatal error occurs, the program\n"
32 "exits with status 1; otherwise, it exits with status 0.\n"
33 "\n"
34 "OPTIONS:\n"
35 " -apar aaa = Use the AFNI dataset 'aaa' as the source of the\n"
36 " transformation; this dataset must be in +acpc\n"
37 " or +tlrc coordinates, and must contain the\n"
38 " attributes WARP_TYPE and WARP_DATA which describe\n"
39 " the forward transformation from +orig coordinates\n"
40 " to the 'aaa' coordinate system.\n"
41 " N.B.: The +orig version of this dataset must also be\n"
42 " readable, since it is also needed when translating\n"
43 " vectors between SureFit and AFNI coordinates.\n"
44 " Only the .HEAD files are actually used.\n"
45 "\n"
46 " -matvec mmm = Read an affine transformation matrix-vector from file\n"
47 " 'mmm', which must be in the format\n"
48 " u11 u12 u13 v1\n"
49 " u21 u22 u23 v2\n"
50 " u31 u32 u33 v3\n"
51 " where each 'uij' and 'vi' is a number. The forward\n"
52 " transformation is defined as\n"
53 " [ xout ] [ u11 u12 u13 ] [ xin ] [ v1 ]\n"
54 " [ yout ] = [ u21 u22 u23 ] [ yin ] + [ v2 ]\n"
55 " [ zout ] [ u31 u32 u33 ] [ zin ] [ v3 ]\n"
56 "\n"
57 " Exactly one of -apar or -matvec must be used to specify the\n"
58 " transformation.\n"
59 "\n"
60 " -forward = -forward means to apply the forward transformation;\n"
61 " *OR* -backward means to apply the backward transformation\n"
62 " -backward * For example, if the transformation is specified by\n"
63 " '-apar fred+tlrc', then the forward transformation\n"
64 " is from +orig to +tlrc coordinates, and the backward\n"
65 " transformation is from +tlrc to +orig coordinates.\n"
66 " * If the transformation is specified by -matvec, then\n"
67 " the matrix-vector read in defines the forward\n"
68 " transform as above, and the backward transformation\n"
69 " is defined as the inverse.\n"
70 " * If neither -forward nor -backward is given, then\n"
71 " -forward is the default.\n"
72 "\n"
73 " -input iii = Read input 3-vectors from file 'iii' (from stdin if\n"
74 " 'iii' is '-' or the -input option is missing). Input\n"
75 " data may be in one of the following ASCII formats:\n"
76 "\n"
77 " * SureFit .coord files:\n"
78 " BeginHeader\n"
79 " lines of text ...\n"
80 " EndHeader\n"
81 " count\n"
82 " int x y z\n"
83 " int x y z\n"
84 " et cetera...\n"
85 " In this case, everything up to and including the\n"
86 " count is simply passed through to the output. Each\n"
87 " (x,y,z) triple is transformed, and output with the\n"
88 " int label that precedes it. Lines that cannot be\n"
89 " scanned as 1 int and 3 floats are treated as comments\n"
90 " and are passed to through to the output unchanged.\n"
91 " N.B.-1: For those using SureFit surfaces created after\n"
92 " the SureFit/Caret merger (post. 2005), you need\n"
93 " to use the flag -new_surefit. Talk to Donna about\n"
94 " this!\n"
95 " N.B.-2: SureFit coordinates are\n"
96 " x = distance Right of Left-most dataset corner\n"
97 " y = distance Anterior to Posterior-most dataset corner\n"
98 " z = distance Superior to Inferior-most dataset corner\n"
99 " For example, if the transformation is specified by\n"
100 " -forward -apar fred+tlrc\n"
101 " then the input (x,y,z) are relative to fred+orig and the\n"
102 " output (x,y,z) are relative to fred+tlrc. If instead\n"
103 " -backward -apar fred+tlrc\n"
104 " is used, then the input (x,y,z) are relative to fred+tlrc\n"
105 " and the output (x,y,z) are relative to fred+orig.\n"
106 " For this to work properly, not only fred+tlrc must be\n"
107 " readable by Vecwarp, but fred+orig must be as well.\n"
108 " If the transformation is specified by -matvec, then\n"
109 " the matrix-vector transformation is applied to the\n"
110 " (x,y,z) vectors directly, with no coordinate shifting.\n"
111 "\n"
112 " * AFNI .1D files with 3 columns\n"
113 " x y z\n"
114 " x y z\n"
115 " et cetera...\n"
116 " In this case, each (x,y,z) triple is transformed and\n"
117 " written to the output. Lines that cannot be scanned\n"
118 " as 3 floats are treated as comments and are passed\n"
119 " through to the output unchanged.\n"
120 " N.B.: AFNI (x,y,z) coordinates are in DICOM order:\n"
121 " -x = Right +x = Left\n"
122 " -y = Anterior +y = Posterior\n"
123 " -z = Inferior +z = Superior\n"
124 "\n"
125 " -output ooo = Write the output to file 'ooo' (to stdout if 'ooo'\n"
126 " is '-', or if the -output option is missing). If the\n"
127 " file already exists, it will not be overwritten unless\n"
128 " the -force option is also used.\n"
129 "\n"
130 " -force = If the output file already exists, -force can be\n"
131 " used to overwrite it. If you want to use -force,\n"
132 " it must come before -output on the command line.\n"
133 "\n"
134 "EXAMPLES:\n"
135 "\n"
136 " Vecwarp -apar fred+tlrc -input fred.orig.coord > fred.tlrc.coord\n"
137 "\n"
138 "This transforms the vectors defined in original coordinates to\n"
139 "Talairach coordinates, using the transformation previously defined\n"
140 "by AFNI markers.\n"
141 "\n"
142 " Vecwarp -apar fred+tlrc -input fred.tlrc.coord -backward > fred.test.coord\n"
143 "\n"
144 "This does the reverse transformation; fred.test.coord should differ from\n"
145 "fred.orig.coord only by roundoff error.\n"
146 "\n"
147 "Author: RWCox - October 2001\n"
148 ) ;
149 PRINT_COMPILE_DATE ; exit(0) ;
150 }
151
152 /*--------------------------------------------------------------------------*/
153
errex(char * str)154 static void errex( char * str )
155 {
156 fprintf(stderr,"** FATAL ERROR: %s\n",str) ; exit(1) ;
157 }
158
159 /*--------------------------------------------------------------------------*/
160
161 #define NBUF 1024
162 #define SUREFIT 33
163 #define AFNI_1D 44
164
main(int argc,char * argv[])165 int main( int argc , char *argv[] )
166 {
167 int iarg=1 ;
168 FILE *fin=stdin , *fout=stdout ;
169 int backward=0 , force=0 ;
170 THD_warp *warp=NULL ;
171 char lbuf[NBUF] , *cpt ;
172 int itype=AFNI_1D , numv=0 , numc=0 ;
173 THD_fvec3 vin , vout ;
174 float xx,yy,zz ;
175 int nn , ii , good=0 ;
176 byte old_surefit_mode = 1;
177 THD_3dim_dataset *aset=NULL , *oset=NULL ; char *oname=NULL ;
178
179 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Syntax() ;
180
181 /*-- process command line arguments --*/
182 old_surefit_mode = 1;
183 while( iarg < argc ){
184
185 /* -input */
186
187 if( strcmp(argv[iarg],"-input") == 0 ){
188 if( ++iarg >= argc )
189 errex("-input: Need argument after -input") ;
190 if( fin != stdin )
191 errex("-input: Can't have two -input options") ;
192 if( strcmp(argv[iarg],"-") != 0 ){
193 fin = fopen( argv[iarg] , "r" ) ;
194 if( fin == NULL )
195 errex("-input: Can't open input file") ;
196 }
197 iarg++ ; continue ;
198 }
199
200 /* -output */
201
202 if( strcmp(argv[iarg],"-output") == 0 ){
203 if( ++iarg >= argc )
204 errex("-output: Need argument after -output") ;
205 if( fout != stdout )
206 errex("-output: Can't have two -output options") ;
207 if( strcmp(argv[iarg],"-") != 0 ){
208 if( !THD_filename_ok(argv[iarg]) )
209 errex("-output: Output filename is illegal") ;
210 if( !force && THD_is_file(argv[iarg]) )
211 errex("-output: Output file already exists") ;
212 fout = fopen( argv[iarg] , "w" ) ;
213 if( fout == NULL )
214 errex("-output: Can't open output file") ;
215 }
216 iarg++ ; continue ;
217 }
218
219 /* -force */
220
221 if( strcmp(argv[iarg],"-force") == 0 ){
222 force = 1 ;
223 iarg++ ; continue ;
224 }
225
226 /* -new_surefit */
227
228 if( strcmp(argv[iarg],"-new_surefit") == 0 ){
229 old_surefit_mode = 0 ;
230 iarg++ ; continue ;
231 }
232
233 /* -forward */
234
235 if( strcmp(argv[iarg],"-forward") == 0 ){
236 backward = 0 ;
237 iarg++ ; continue ;
238 }
239
240 /* -backward */
241
242 if( strcmp(argv[iarg],"-backward") == 0 ){
243 backward = 1 ;
244 iarg++ ; continue ;
245 }
246
247 /* -apar */
248
249 if( strcmp(argv[iarg],"-apar") == 0 ){
250 if( ++iarg >= argc )
251 errex("-apar: Need argument after -apar") ;
252 if( warp != NULL )
253 errex("-apar: Can't specify transformation twice") ;
254
255 /* open dataset with warp */
256
257 aset = THD_open_dataset( argv[iarg] ) ;
258 if( !ISVALID_DSET(aset) ){
259 sprintf(lbuf,"-apar: Can't open dataset %s\n",argv[iarg]) ;
260 errex(lbuf) ;
261 }
262 if( aset->warp == NULL ){
263 sprintf(lbuf,"-apar: Dataset %s does not contain warp",argv[iarg]) ;
264 errex(lbuf) ;
265 }
266 if( aset->view_type == VIEW_ORIGINAL_TYPE ){
267 sprintf(lbuf,"-apar: Dataset %s is in the +orig view",argv[iarg]) ;
268 errex(lbuf) ;
269 }
270
271 /* open +orig version of this dataset */
272
273 sprintf(lbuf,"%s%s+orig.HEAD", aset->dblk->diskptr->directory_name,
274 aset->dblk->diskptr->prefix );
275
276 oname = strdup(lbuf) ;
277 oset = THD_open_dataset(lbuf) ;
278 #if 0
279 if( !ISVALID_DSET(oset) ){
280 char str[NBUF] ;
281 sprintf(str,"-apar: Can't open dataset %s",lbuf) ;
282 errex(str) ;
283 }
284 #endif
285
286 warp = aset->warp ;
287 iarg++ ; continue ;
288 }
289
290 /* -matvec */
291
292 if( strcmp(argv[iarg],"-matvec") == 0 ){
293 THD_dvecmat dvm ;
294 THD_linear_mapping lmap ;
295
296 ZZME(lmap) ;
297
298 if( ++iarg >= argc )
299 errex("-matvec: Need argument after -matvec") ;
300 if( warp != NULL )
301 errex("-matvec: Can't specify transformation twice") ;
302 dvm = THD_read_dvecmat( argv[iarg] , 0 ) ;
303 if( SIZE_DMAT(dvm.mm) == 0.0 )
304 errex("-matvec: Can't read transformation") ;
305
306 /* manufacture an AFNI linear map */
307
308 lmap.type = MAPPING_LINEAR_TYPE ;
309 lmap.mfor.mat[0][0] = dvm.mm.mat[0][0] ; /* copy matrix in */
310 lmap.mfor.mat[0][1] = dvm.mm.mat[0][1] ;
311 lmap.mfor.mat[0][2] = dvm.mm.mat[0][2] ;
312 lmap.mfor.mat[1][0] = dvm.mm.mat[1][0] ;
313 lmap.mfor.mat[1][1] = dvm.mm.mat[1][1] ;
314 lmap.mfor.mat[1][2] = dvm.mm.mat[1][2] ;
315 lmap.mfor.mat[2][0] = dvm.mm.mat[2][0] ;
316 lmap.mfor.mat[2][1] = dvm.mm.mat[2][1] ;
317 lmap.mfor.mat[2][2] = dvm.mm.mat[2][2] ;
318 lmap.bvec.xyz[0] = -dvm.vv.xyz[0] ; /* copy vector in */
319 lmap.bvec.xyz[1] = -dvm.vv.xyz[1] ;
320 lmap.bvec.xyz[2] = -dvm.vv.xyz[2] ;
321 LOAD_INVERSE_LMAP(lmap) ; /* make inverse map */
322
323 /* manufacture an AFNI warp */
324
325 warp = (THD_warp *) calloc( sizeof(THD_warp) , 1 ) ;
326 warp->type = WARP_AFFINE_TYPE ;
327 warp->rig_bod.warp = lmap ;
328
329 /* corrected the output location for vv.xyz 17 Aug 2006 [rickr] */
330 fprintf(stderr,"++ Using matrix-vector transformation below:\n"
331 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n"
332 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n"
333 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n" ,
334 dvm.mm.mat[0][0] , dvm.mm.mat[0][1] , dvm.mm.mat[0][2] , dvm.vv.xyz[0] ,
335 dvm.mm.mat[1][0] , dvm.mm.mat[1][1] , dvm.mm.mat[1][2] , dvm.vv.xyz[1] ,
336 dvm.mm.mat[2][0] , dvm.mm.mat[2][1] , dvm.mm.mat[2][2] , dvm.vv.xyz[2]) ;
337
338 iarg++ ; continue ;
339 }
340
341 /* ERROR */
342
343 { char *str = AFMALL(char, strlen(argv[iarg])+32) ;
344 sprintf(str,"Unknown option: %s",argv[iarg]) ;
345 errex(str) ;
346 }
347
348 } /* end of loop over command line args */
349
350 /*-- check if we are ready --*/
351
352 if( warp == NULL )
353 errex("Transformation not specified") ;
354
355 /*-- Read 1st line of input file to determine what to do with it --*/
356
357 cpt = afni_fgets( lbuf , NBUF , fin ) ;
358 if( cpt == NULL )
359 errex("Couldn't read 1st line from input file") ;
360
361 /*-- check if this is a SureFit file --*/
362
363 if( strstr(cpt,"BeginHeader") != NULL ) itype = SUREFIT ;
364 else itype = AFNI_1D ;
365
366 if( itype == SUREFIT && oset == NULL ){
367 ERROR_exit(
368 "SureFit transformation requires that the +orig dataset also exists:\n"
369 " %s" ,
370 oname ) ;
371 }
372
373 /*-- if SureFit, echo all header stuff to the output --*/
374
375 if( itype == SUREFIT ){
376 int numh=1 ;
377 fprintf(fout,"%s",lbuf) ; /* echo 1st line */
378 while(1){ /* read lines until EndHeader is found */
379
380 cpt = afni_fgets( lbuf , NBUF , fin ) ;
381 if( cpt == NULL )
382 errex("Input ended before EndHeader was found") ;
383
384 fprintf(fout,"%s",lbuf) ; /* echo line to output */
385 numh++ ;
386
387 if( strstr(lbuf,"EndHeader") != NULL ){ /* do next line, too */
388 cpt = afni_fgets( lbuf , NBUF , fin ) ;
389 if( cpt == NULL )
390 errex("Input ended just after EndHeader") ;
391 fprintf(fout,"%s",lbuf) ;
392 fprintf(stderr,"++ Wrote %d SureFit header lines to output\n",numh);
393 break ; /* end of loop */
394 }
395 }
396
397 cpt = afni_fgets( lbuf , NBUF , fin ) ; /* get next line, with 1st vector */
398 if( cpt == NULL )
399 errex("Input ended just after Node count") ;
400
401 } /* end of SureFit header echo */
402
403 /*-- From this point on, process each line with 1 vector. --*/
404 /*-- At the start of the loop, 1 line's data is in lbuf. --*/
405
406 do{
407
408 switch( itype ){
409 case SUREFIT:
410 ii = sscanf(lbuf,"%d%f%f%f",&nn,&xx,&yy,&zz) ;
411 good = (ii == 4) ;
412 break ;
413
414 case AFNI_1D:
415 ii = sscanf(lbuf,"%f%f%f",&xx,&yy,&zz) ;
416 good = (ii == 3) ;
417 break ;
418 }
419
420 if( !good ){ /* just echo line to output */
421
422 fprintf(fout,"%s",lbuf) ;
423 numc++ ;
424
425 } else { /* transform a vector!!! */
426
427 LOAD_FVEC3(vin,xx,yy,zz) ;
428 #ifdef DEBUG
429 fprintf(stderr,"\n") ;
430 DUMP_FVEC3("vin ",vin) ;
431 #endif
432
433 if( itype == SUREFIT && old_surefit_mode && aset != NULL ){
434 if( backward ) vin = THD_surefit_to_dicomm( aset , vin ) ;
435 else vin = THD_surefit_to_dicomm( oset , vin ) ;
436 #ifdef DEBUG
437 DUMP_FVEC3("surefit_to_dicomm",vin) ;
438 #endif
439 }
440
441 if( backward ) vout = myAFNI_backward_warp_vector( warp , vin ) ;
442 else vout = myAFNI_forward_warp_vector ( warp , vin ) ;
443
444 #ifdef DEBUG
445 DUMP_FVEC3("vout ",vout) ;
446 #endif
447
448 if( itype == SUREFIT && old_surefit_mode && aset != NULL ){
449 if( backward ) vout = THD_dicomm_to_surefit( oset , vout ) ;
450 else vout = THD_dicomm_to_surefit( aset , vout ) ;
451 #ifdef DEBUG
452 DUMP_FVEC3("dicomm_to_surefit",vout) ;
453 #endif
454 }
455
456 if( itype == SUREFIT ) fprintf(fout,"%d ",nn) ;
457 fprintf(fout,"%f %f %f\n",vout.xyz[0],vout.xyz[1],vout.xyz[2]) ;
458 numv++ ;
459 }
460
461 cpt = afni_fgets( lbuf , NBUF , fin ) ; /* get next line */
462
463 } while( cpt != NULL ) ; /* loop until no data can be read */
464
465 if( numc > 0 )
466 fprintf(stderr,"++ Wrote %d vectors; %d comments\n",numv,numc) ;
467 else
468 fprintf(stderr,"++ Wrote %d vectors\n",numv) ;
469
470 exit(0) ;
471 }
472
473 /*--------------------------------------------------------------------------
474 The following functions are adapted from afni.c
475 ----------------------------------------------------------------------------*/
476
477 /*------------------------------------------------------------------------
478 Forward transform a vector following a warp
479 --------------------------------------------------------------------------*/
480
myAFNI_forward_warp_vector(THD_warp * warp,THD_fvec3 old_fv)481 THD_fvec3 myAFNI_forward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
482 {
483 THD_fvec3 new_fv ;
484
485 if( warp == NULL ) return old_fv ;
486
487 switch( warp->type ){
488
489 default: new_fv = old_fv ; break ;
490
491 case WARP_TALAIRACH_12_TYPE:{
492 THD_linear_mapping map ;
493 int iw ;
494
495 /* forward transform each possible case,
496 and test if result is in bot..top of defined map */
497
498 for( iw=0 ; iw < 12 ; iw++ ){
499 map = warp->tal_12.warp[iw] ;
500 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
501
502 if( new_fv.xyz[0] >= map.bot.xyz[0] &&
503 new_fv.xyz[1] >= map.bot.xyz[1] &&
504 new_fv.xyz[2] >= map.bot.xyz[2] &&
505 new_fv.xyz[0] <= map.top.xyz[0] &&
506 new_fv.xyz[1] <= map.top.xyz[1] &&
507 new_fv.xyz[2] <= map.top.xyz[2] ) break ; /* leave loop */
508 }
509 }
510 break ;
511
512 case WARP_AFFINE_TYPE:{
513 THD_linear_mapping map = warp->rig_bod.warp ;
514 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
515 }
516 break ;
517
518 }
519 return new_fv ;
520 }
521
522 /*------------------------------------------------------------------------
523 Backward transform a vector following a warp
524 --------------------------------------------------------------------------*/
525
myAFNI_backward_warp_vector(THD_warp * warp,THD_fvec3 old_fv)526 THD_fvec3 myAFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
527 {
528 THD_fvec3 new_fv ;
529
530 if( warp == NULL ) return old_fv ;
531
532 switch( warp->type ){
533
534 default: new_fv = old_fv ; break ;
535
536 case WARP_TALAIRACH_12_TYPE:{
537 THD_linear_mapping map ;
538 int iw ;
539
540 /* test if input is in bot..top of each defined map */
541
542 for( iw=0 ; iw < 12 ; iw++ ){
543 map = warp->tal_12.warp[iw] ;
544
545 if( old_fv.xyz[0] >= map.bot.xyz[0] &&
546 old_fv.xyz[1] >= map.bot.xyz[1] &&
547 old_fv.xyz[2] >= map.bot.xyz[2] &&
548 old_fv.xyz[0] <= map.top.xyz[0] &&
549 old_fv.xyz[1] <= map.top.xyz[1] &&
550 old_fv.xyz[2] <= map.top.xyz[2] ) break ; /* leave loop */
551 }
552 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
553 }
554 break ;
555
556 case WARP_AFFINE_TYPE:{
557 THD_linear_mapping map = warp->rig_bod.warp ;
558 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
559 }
560 break ;
561
562 }
563 return new_fv ;
564 }
565
566 #if 0
567 /*--------------------------------------------------------------------------
568 The following routines are used to convert DICOM order coordinates
569 (used in AFNI) to SureFit order coordinates
570 ----------------------------------------------------------------------------*/
571
572 static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *dset , THD_fvec3 fv )
573 {
574 float xx,yy,zz , xbase,ybase,zbase ;
575 THD_fvec3 vout ;
576 #ifdef DEBUG
577 static int first=1 ;
578 #endif
579
580 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ; /* xyz now LPI */
581
582 if( dset != NULL ){
583 THD_fvec3 v1 , v2 ;
584 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
585 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
586 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
587 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
588 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
589 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
590 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ; /* Left-most */
591 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ; /* Posterior */
592 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ; /* Inferior */
593 } else {
594 xbase = ybase = zbase = 0.0 ;
595 }
596 #ifdef DEBUG
597 if(first){fprintf(stderr,"d2s base=%f %f %f\n",xbase,ybase,zbase);first=0;}
598 #endif
599
600 vout.xyz[0] = xx - xbase ;
601 vout.xyz[1] = yy - ybase ;
602 vout.xyz[2] = zz - zbase ; return vout ;
603 }
604
605 static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 fv )
606 {
607 float xx,yy,zz , xbase,ybase,zbase ;
608 THD_fvec3 vout ;
609 #ifdef DEBUG
610 static int first=1 ;
611 #endif
612
613 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ; /* xyz now RAI */
614
615 if( dset != NULL ){
616 THD_fvec3 v1 , v2 ;
617 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
618 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
619 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
620 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
621 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
622 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
623 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
624 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
625 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
626 } else {
627 xbase = ybase = zbase = 0.0 ;
628 }
629 #ifdef DEBUG
630 if(first){fprintf(stderr,"s2d base=%f %f %f\n",xbase,ybase,zbase);first=0;}
631 #endif
632
633 vout.xyz[0] = xx - xbase ;
634 vout.xyz[1] = yy - ybase ;
635 vout.xyz[2] = zz + zbase ; return vout ;
636 }
637 #endif
638