00001 #include "mrilib.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 static THD_fvec3 AFNI_forward_warp_vector ( THD_warp * , THD_fvec3 ) ;
00014 static THD_fvec3 AFNI_backward_warp_vector( THD_warp * , THD_fvec3 ) ;
00015
00016 #if 0
00017 static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *, THD_fvec3 ) ;
00018 static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *, THD_fvec3 ) ;
00019 #endif
00020
00021 #undef DEBUG
00022
00023
00024
00025 static void Syntax(void)
00026 {
00027 printf(
00028 "Usage: Vecwarp [options]\n"
00029 "Transforms (warps) a list of 3-vectors into another list of 3-vectors\n"
00030 "according to the options. Error messages, warnings, and informational\n"
00031 "messages are written to stderr. If a fatal error occurs, the program\n"
00032 "exits with status 1; otherwise, it exits with status 0.\n"
00033 "\n"
00034 "OPTIONS:\n"
00035 " -apar aaa = Use the AFNI dataset 'aaa' as the source of the\n"
00036 " transformation; this dataset must be in +acpc\n"
00037 " or +tlrc coordinates, and must contain the\n"
00038 " attributes WARP_TYPE and WARP_DATA which describe\n"
00039 " the forward transformation from +orig coordinates\n"
00040 " to the 'aaa' coordinate system.\n"
00041 " N.B.: The +orig version of this dataset must also be\n"
00042 " readable, since it is also needed when translating\n"
00043 " vectors between SureFit and AFNI coordinates.\n"
00044 " Only the .HEAD files are actually used.\n"
00045 "\n"
00046 " -matvec mmm = Read an affine transformation matrix-vector from file\n"
00047 " 'mmm', which must be in the format\n"
00048 " u11 u12 u13 v1\n"
00049 " u21 u22 u23 v2\n"
00050 " u31 u32 u33 v3\n"
00051 " where each 'uij' and 'vi' is a number. The forward\n"
00052 " transformation is defined as\n"
00053 " [ xout ] [ u11 u12 u13 ] [ xin ] [ v1 ]\n"
00054 " [ yout ] = [ u21 u22 u23 ] [ yin ] + [ v2 ]\n"
00055 " [ zout ] [ u31 u32 u33 ] [ zin ] [ v3 ]\n"
00056 "\n"
00057 " Exactly one of -apar or -matvec must be used to specify the\n"
00058 " transformation.\n"
00059 "\n"
00060 " -forward = -forward means to apply the forward transformation;\n"
00061 " *OR* -backward means to apply the backward transformation\n"
00062 " -backward * For example, if the transformation is specified by\n"
00063 " '-apar fred+tlrc', then the forward transformation\n"
00064 " is from +orig to +tlrc coordinates, and the backward\n"
00065 " transformation is from +tlrc to +orig coordinates.\n"
00066 " * If the transformation is specified by -matvec, then\n"
00067 " the matrix-vector read in defines the forward\n"
00068 " transform as above, and the backward transformation\n"
00069 " is defined as the inverse.\n"
00070 " * If neither -forward nor -backward is given, then\n"
00071 " -forward is the default.\n"
00072 "\n"
00073 " -input iii = Read input 3-vectors from file 'iii' (from stdin if\n"
00074 " 'iii' is '-' or the -input option is missing). Input\n"
00075 " data may be in one of the following ASCII formats:\n"
00076 "\n"
00077 " * SureFit .coord files:\n"
00078 " BeginHeader\n"
00079 " lines of text ...\n"
00080 " EndHeader\n"
00081 " count\n"
00082 " int x y z\n"
00083 " int x y z\n"
00084 " et cetera...\n"
00085 " In this case, everything up to and including the\n"
00086 " count is simply passed through to the output. Each\n"
00087 " (x,y,z) triple is transformed, and output with the\n"
00088 " int label that precedes it. Lines that cannot be\n"
00089 " scanned as 1 int and 3 floats are treated as comments\n"
00090 " and are passed to through to the output unchanged.\n"
00091 " N.B.: SureFit coordinates are\n"
00092 " x = distance Right of Left-most dataset corner\n"
00093 " y = distance Anterior to Posterior-most dataset corner\n"
00094 " z = distance Superior to Inferior-most dataset corner\n"
00095 " For example, if the transformation is specified by\n"
00096 " -forward -apar fred+tlrc\n"
00097 " then the input (x,y,z) are relative to fred+orig and the\n"
00098 " output (x,y,z) are relative to fred+tlrc. If instead\n"
00099 " -backward -apar fred+tlrc\n"
00100 " is used, then the input (x,y,z) are relative to fred+tlrc\n"
00101 " and the output (x,y,z) are relative to fred+orig.\n"
00102 " For this to work properly, not only fred+tlrc must be\n"
00103 " readable by Vecwarp, but fred+orig must be as well.\n"
00104 " If the transformation is specified by -matvec, then\n"
00105 " the matrix-vector transformation is applied to the\n"
00106 " (x,y,z) vectors directly, with no coordinate shifting.\n"
00107 "\n"
00108 " * AFNI .1D files with 3 columns\n"
00109 " x y z\n"
00110 " x y z\n"
00111 " et cetera...\n"
00112 " In this case, each (x,y,z) triple is transformed and\n"
00113 " written to the output. Lines that cannot be scanned\n"
00114 " as 3 floats are treated as comments and are passed\n"
00115 " through to the output unchanged.\n"
00116 " N.B.: AFNI (x,y,z) coordinates are in DICOM order:\n"
00117 " -x = Right +x = Left\n"
00118 " -y = Anterior +y = Posterior\n"
00119 " -z = Inferior +z = Superior\n"
00120 "\n"
00121 " -output ooo = Write the output to file 'ooo' (to stdout if 'ooo'\n"
00122 " is '-', or if the -output option is missing). If the\n"
00123 " file already exists, it will not be overwritten unless\n"
00124 " the -force option is also used.\n"
00125 "\n"
00126 " -force = If the output file already exists, -force can be\n"
00127 " used to overwrite it. If you want to use -force,\n"
00128 " it must come before -output on the command line.\n"
00129 "\n"
00130 "EXAMPLES:\n"
00131 "\n"
00132 " Vecwarp -apar fred+tlrc -input fred.orig.coord > fred.tlrc.coord\n"
00133 "\n"
00134 "This transforms the vectors defined in original coordinates to\n"
00135 "Talairach coordinates, using the transformation previously defined\n"
00136 "by AFNI markers.\n"
00137 "\n"
00138 " Vecwarp -apar fred+tlrc -input fred.tlrc.coord -backward > fred.test.coord\n"
00139 "\n"
00140 "This does the reverse transformation; fred.test.coord should differ from\n"
00141 "fred.orig.coord only by roundoff error.\n"
00142 "\n"
00143 "Author: RWCox - October 2001\n"
00144 ) ;
00145 exit(0) ;
00146 }
00147
00148
00149
00150 static void errex( char * str )
00151 {
00152 fprintf(stderr,"** FATAL ERROR: %s\n",str) ; exit(1) ;
00153 }
00154
00155
00156
00157 #define NBUF 1024
00158 #define isnumeric(c) (isdigit(c) || c == '-' || c == '+' || c == '.')
00159
00160 #define SUREFIT 33
00161 #define AFNI_1D 44
00162
00163 int main( int argc , char *argv[] )
00164 {
00165 int iarg=1 ;
00166 FILE *fin=stdin , *fout=stdout ;
00167 int backward=0 , force=0 ;
00168 THD_warp *warp=NULL ;
00169 char lbuf[NBUF] , *cpt ;
00170 int itype=AFNI_1D , numv=0 , numc=0 ;
00171 THD_fvec3 vin , vout ;
00172 float xx,yy,zz ;
00173 int nn , ii , good=0 ;
00174 THD_3dim_dataset *aset=NULL , *oset=NULL ;
00175
00176 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Syntax() ;
00177
00178
00179
00180 while( iarg < argc ){
00181
00182
00183
00184 if( strcmp(argv[iarg],"-input") == 0 ){
00185 if( ++iarg >= argc )
00186 errex("-input: Need argument after -input") ;
00187 if( fin != stdin )
00188 errex("-input: Can't have two -input options") ;
00189 if( strcmp(argv[iarg],"-") != 0 ){
00190 fin = fopen( argv[iarg] , "r" ) ;
00191 if( fin == NULL )
00192 errex("-input: Can't open input file") ;
00193 }
00194 iarg++ ; continue ;
00195 }
00196
00197
00198
00199 if( strcmp(argv[iarg],"-output") == 0 ){
00200 if( ++iarg >= argc )
00201 errex("-output: Need argument after -output") ;
00202 if( fout != stdout )
00203 errex("-output: Can't have two -output options") ;
00204 if( strcmp(argv[iarg],"-") != 0 ){
00205 if( !THD_filename_ok(argv[iarg]) )
00206 errex("-output: Output filename is illegal") ;
00207 if( !force && THD_is_file(argv[iarg]) )
00208 errex("-output: Output file already exists") ;
00209 fout = fopen( argv[iarg] , "w" ) ;
00210 if( fout == NULL )
00211 errex("-output: Can't open output file") ;
00212 }
00213 iarg++ ; continue ;
00214 }
00215
00216
00217
00218 if( strcmp(argv[iarg],"-force") == 0 ){
00219 force = 1 ;
00220 iarg++ ; continue ;
00221 }
00222
00223
00224
00225 if( strcmp(argv[iarg],"-forward") == 0 ){
00226 backward = 0 ;
00227 iarg++ ; continue ;
00228 }
00229
00230
00231
00232 if( strcmp(argv[iarg],"-backward") == 0 ){
00233 backward = 1 ;
00234 iarg++ ; continue ;
00235 }
00236
00237
00238
00239 if( strcmp(argv[iarg],"-apar") == 0 ){
00240 if( ++iarg >= argc )
00241 errex("-apar: Need argument after -apar") ;
00242 if( warp != NULL )
00243 errex("-apar: Can't specify transformation twice") ;
00244
00245
00246
00247 aset = THD_open_dataset( argv[iarg] ) ;
00248 if( !ISVALID_DSET(aset) ){
00249 sprintf(lbuf,"-apar: Can't open dataset %s\n",argv[iarg]) ;
00250 errex(lbuf) ;
00251 }
00252 if( aset->warp == NULL ){
00253 sprintf(lbuf,"-apar: Dataset %s does not contain warp",argv[iarg]) ;
00254 errex(lbuf) ;
00255 }
00256 if( aset->view_type == VIEW_ORIGINAL_TYPE ){
00257 sprintf(lbuf,"-apar: Dataset %s is in the +orig view",argv[iarg]) ;
00258 errex(lbuf) ;
00259 }
00260
00261
00262
00263 sprintf(lbuf,"%s%s+orig.HEAD", aset->dblk->diskptr->directory_name,
00264 aset->dblk->diskptr->prefix );
00265
00266 oset = THD_open_dataset(lbuf) ;
00267 if( !ISVALID_DSET(oset) ){
00268 char str[NBUF] ;
00269 sprintf(str,"-apar: Can't open dataset %s",lbuf) ;
00270 errex(str) ;
00271 }
00272
00273 warp = aset->warp ;
00274 iarg++ ; continue ;
00275 }
00276
00277
00278
00279 if( strcmp(argv[iarg],"-matvec") == 0 ){
00280 THD_dvecmat dvm ;
00281 THD_linear_mapping lmap ;
00282
00283 if( ++iarg >= argc )
00284 errex("-matvec: Need argument after -matvec") ;
00285 if( warp != NULL )
00286 errex("-matvec: Can't specify transformation twice") ;
00287 dvm = THD_read_dvecmat( argv[iarg] , 0 ) ;
00288 if( SIZE_DMAT(dvm.mm) == 0.0 )
00289 errex("-matvec: Can't read transformation") ;
00290
00291
00292
00293 lmap.type = MAPPING_LINEAR_TYPE ;
00294 lmap.mfor.mat[0][0] = dvm.mm.mat[0][0] ;
00295 lmap.mfor.mat[0][1] = dvm.mm.mat[0][1] ;
00296 lmap.mfor.mat[0][2] = dvm.mm.mat[0][2] ;
00297 lmap.mfor.mat[1][0] = dvm.mm.mat[1][0] ;
00298 lmap.mfor.mat[1][1] = dvm.mm.mat[1][1] ;
00299 lmap.mfor.mat[1][2] = dvm.mm.mat[1][2] ;
00300 lmap.mfor.mat[2][0] = dvm.mm.mat[2][0] ;
00301 lmap.mfor.mat[2][1] = dvm.mm.mat[2][1] ;
00302 lmap.mfor.mat[2][2] = dvm.mm.mat[2][2] ;
00303 lmap.bvec.xyz[0] = -dvm.vv.xyz[0] ;
00304 lmap.bvec.xyz[1] = -dvm.vv.xyz[1] ;
00305 lmap.bvec.xyz[2] = -dvm.vv.xyz[2] ;
00306 LOAD_INVERSE_LMAP(lmap) ;
00307
00308
00309
00310 warp = (THD_warp *) calloc( sizeof(THD_warp) , 1 ) ;
00311 warp->type = WARP_AFFINE_TYPE ;
00312 warp->rig_bod.warp = lmap ;
00313
00314 fprintf(stderr,"++ Using matrix-vector transformation below:\n"
00315 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n"
00316 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n"
00317 " [ %9.5f %9.5f %9.5f ] [ %9.5f ]\n" ,
00318 dvm.mm.mat[0][0] , dvm.mm.mat[0][1] , dvm.mm.mat[0][2] ,
00319 dvm.mm.mat[1][0] , dvm.mm.mat[1][1] , dvm.mm.mat[1][2] ,
00320 dvm.mm.mat[2][0] , dvm.mm.mat[2][1] , dvm.mm.mat[2][2] ,
00321 dvm.vv.xyz[0] , dvm.vv.xyz[1] , dvm.vv.xyz[2] ) ;
00322
00323 iarg++ ; continue ;
00324 }
00325
00326
00327
00328 { char *str = AFMALL(char, strlen(argv[iarg])+32) ;
00329 sprintf(str,"Unknown option: %s",argv[iarg]) ;
00330 errex(str) ;
00331 }
00332
00333 }
00334
00335
00336
00337 if( warp == NULL )
00338 errex("Transformation not specified") ;
00339
00340
00341
00342 cpt = fgets( lbuf , NBUF , fin ) ;
00343 if( cpt == NULL )
00344 errex("Couldn't read 1st line from input file") ;
00345
00346
00347
00348 if( strstr(cpt,"BeginHeader") != NULL ) itype = SUREFIT ;
00349 else itype = AFNI_1D ;
00350
00351
00352
00353 if( itype == SUREFIT ){
00354 int numh=1 ;
00355 fprintf(fout,"%s",lbuf) ;
00356 while(1){
00357
00358 cpt = fgets( lbuf , NBUF , fin ) ;
00359 if( cpt == NULL )
00360 errex("Input ended before EndHeader was found") ;
00361
00362 fprintf(fout,"%s",lbuf) ;
00363 numh++ ;
00364
00365 if( strstr(lbuf,"EndHeader") != NULL ){
00366 cpt = fgets( lbuf , NBUF , fin ) ;
00367 if( cpt == NULL )
00368 errex("Input ended just after EndHeader") ;
00369 fprintf(fout,"%s",lbuf) ;
00370 fprintf(stderr,"++ Wrote %d SureFit header lines to output\n",numh);
00371 break ;
00372 }
00373 }
00374
00375 cpt = fgets( lbuf , NBUF , fin ) ;
00376 if( cpt == NULL )
00377 errex("Input ended just after Node count") ;
00378
00379 }
00380
00381
00382
00383
00384 do{
00385
00386 switch( itype ){
00387 case SUREFIT:
00388 ii = sscanf(lbuf,"%d%f%f%f",&nn,&xx,&yy,&zz) ;
00389 good = (ii == 4) ;
00390 break ;
00391
00392 case AFNI_1D:
00393 ii = sscanf(lbuf,"%f%f%f",&xx,&yy,&zz) ;
00394 good = (ii == 3) ;
00395 break ;
00396 }
00397
00398 if( !good ){
00399
00400 fprintf(fout,"%s",lbuf) ;
00401 numc++ ;
00402
00403 } else {
00404
00405 LOAD_FVEC3(vin,xx,yy,zz) ;
00406 #ifdef DEBUG
00407 fprintf(stderr,"\n") ;
00408 DUMP_FVEC3("vin ",vin) ;
00409 #endif
00410
00411 if( itype == SUREFIT && aset != NULL ){
00412 if( backward ) vin = THD_surefit_to_dicomm( aset , vin ) ;
00413 else vin = THD_surefit_to_dicomm( oset , vin ) ;
00414 #ifdef DEBUG
00415 DUMP_FVEC3("surefit_to_dicomm",vin) ;
00416 #endif
00417 }
00418
00419 if( backward ) vout = AFNI_backward_warp_vector( warp , vin ) ;
00420 else vout = AFNI_forward_warp_vector ( warp , vin ) ;
00421
00422 #ifdef DEBUG
00423 DUMP_FVEC3("vout ",vout) ;
00424 #endif
00425
00426 if( itype == SUREFIT && aset != NULL ){
00427 if( backward ) vout = THD_dicomm_to_surefit( oset , vout ) ;
00428 else vout = THD_dicomm_to_surefit( aset , vout ) ;
00429 #ifdef DEBUG
00430 DUMP_FVEC3("dicomm_to_surefit",vout) ;
00431 #endif
00432 }
00433
00434 if( itype == SUREFIT ) fprintf(fout,"%d ",nn) ;
00435 fprintf(fout,"%f %f %f\n",vout.xyz[0],vout.xyz[1],vout.xyz[2]) ;
00436 numv++ ;
00437 }
00438
00439 cpt = fgets( lbuf , NBUF , fin ) ;
00440
00441 } while( cpt != NULL ) ;
00442
00443 if( numc > 0 )
00444 fprintf(stderr,"++ Wrote %d vectors; %d comments\n",numv,numc) ;
00445 else
00446 fprintf(stderr,"++ Wrote %d vectors\n",numv) ;
00447
00448 exit(0) ;
00449 }
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 static THD_fvec3 AFNI_forward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
00460 {
00461 THD_fvec3 new_fv ;
00462
00463 if( warp == NULL ) return old_fv ;
00464
00465 switch( warp->type ){
00466
00467 default: new_fv = old_fv ; break ;
00468
00469 case WARP_TALAIRACH_12_TYPE:{
00470 THD_linear_mapping map ;
00471 int iw ;
00472
00473
00474
00475
00476 for( iw=0 ; iw < 12 ; iw++ ){
00477 map = warp->tal_12.warp[iw] ;
00478 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00479
00480 if( new_fv.xyz[0] >= map.bot.xyz[0] &&
00481 new_fv.xyz[1] >= map.bot.xyz[1] &&
00482 new_fv.xyz[2] >= map.bot.xyz[2] &&
00483 new_fv.xyz[0] <= map.top.xyz[0] &&
00484 new_fv.xyz[1] <= map.top.xyz[1] &&
00485 new_fv.xyz[2] <= map.top.xyz[2] ) break ;
00486 }
00487 }
00488 break ;
00489
00490 case WARP_AFFINE_TYPE:{
00491 THD_linear_mapping map = warp->rig_bod.warp ;
00492 new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00493 }
00494 break ;
00495
00496 }
00497 return new_fv ;
00498 }
00499
00500
00501
00502
00503
00504 static THD_fvec3 AFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
00505 {
00506 THD_fvec3 new_fv ;
00507
00508 if( warp == NULL ) return old_fv ;
00509
00510 switch( warp->type ){
00511
00512 default: new_fv = old_fv ; break ;
00513
00514 case WARP_TALAIRACH_12_TYPE:{
00515 THD_linear_mapping map ;
00516 int iw ;
00517
00518
00519
00520 for( iw=0 ; iw < 12 ; iw++ ){
00521 map = warp->tal_12.warp[iw] ;
00522
00523 if( old_fv.xyz[0] >= map.bot.xyz[0] &&
00524 old_fv.xyz[1] >= map.bot.xyz[1] &&
00525 old_fv.xyz[2] >= map.bot.xyz[2] &&
00526 old_fv.xyz[0] <= map.top.xyz[0] &&
00527 old_fv.xyz[1] <= map.top.xyz[1] &&
00528 old_fv.xyz[2] <= map.top.xyz[2] ) break ;
00529 }
00530 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00531 }
00532 break ;
00533
00534 case WARP_AFFINE_TYPE:{
00535 THD_linear_mapping map = warp->rig_bod.warp ;
00536 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00537 }
00538 break ;
00539
00540 }
00541 return new_fv ;
00542 }
00543
00544 #if 0
00545
00546
00547
00548
00549
00550 static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *dset , THD_fvec3 fv )
00551 {
00552 float xx,yy,zz , xbase,ybase,zbase ;
00553 THD_fvec3 vout ;
00554 #ifdef DEBUG
00555 static int first=1 ;
00556 #endif
00557
00558 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;
00559
00560 if( dset != NULL ){
00561 THD_fvec3 v1 , v2 ;
00562 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
00563 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
00564 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
00565 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
00566 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
00567 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
00568 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
00569 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
00570 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
00571 } else {
00572 xbase = ybase = zbase = 0.0 ;
00573 }
00574 #ifdef DEBUG
00575 if(first){fprintf(stderr,"d2s base=%f %f %f\n",xbase,ybase,zbase);first=0;}
00576 #endif
00577
00578 vout.xyz[0] = xx - xbase ;
00579 vout.xyz[1] = yy - ybase ;
00580 vout.xyz[2] = zz - zbase ; return vout ;
00581 }
00582
00583 static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 fv )
00584 {
00585 float xx,yy,zz , xbase,ybase,zbase ;
00586 THD_fvec3 vout ;
00587 #ifdef DEBUG
00588 static int first=1 ;
00589 #endif
00590
00591 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;
00592
00593 if( dset != NULL ){
00594 THD_fvec3 v1 , v2 ;
00595 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
00596 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
00597 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
00598 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
00599 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
00600 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
00601 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
00602 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
00603 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
00604 } else {
00605 xbase = ybase = zbase = 0.0 ;
00606 }
00607 #ifdef DEBUG
00608 if(first){fprintf(stderr,"s2d base=%f %f %f\n",xbase,ybase,zbase);first=0;}
00609 #endif
00610
00611 vout.xyz[0] = xx - xbase ;
00612 vout.xyz[1] = yy - ybase ;
00613 vout.xyz[2] = zz + zbase ; return vout ;
00614 }
00615 #endif
