#include "mrilib.h"

/*--------------------------------------------------------------------------
  Program to transform 3-vectors based on warps stored in AFNI .HEAD files.
  Transformations are stored in +view.HEAD files (view=acpc or tlrc) between
  the +orig and +view coordinate systems.

  24 Oct 2001: creation by RWCox
----------------------------------------------------------------------------*/

/*** Some prototypes ***/

THD_fvec3 myAFNI_forward_warp_vector ( THD_warp * , THD_fvec3 ) ;
THD_fvec3 myAFNI_backward_warp_vector( THD_warp * , THD_fvec3 ) ;

#if 0
static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *, THD_fvec3 ) ;
static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *, THD_fvec3 ) ;
#endif

#undef DEBUG

/*--------------------------------------------------------------------------*/

static void Syntax(void)
{
   printf(
    "Usage: Vecwarp [options]\n"
    "Transforms (warps) a list of 3-vectors into another list of 3-vectors\n"
    "according to the options.  Error messages, warnings, and informational\n"
    "messages are written to stderr.  If a fatal error occurs, the program\n"
    "exits with status 1; otherwise, it exits with status 0.\n"
    "\n"
    "OPTIONS:\n"
    " -apar aaa   = Use the AFNI dataset 'aaa' as the source of the\n"
    "               transformation; this dataset must be in +acpc\n"
    "               or +tlrc coordinates, and must contain the\n"
    "               attributes WARP_TYPE and WARP_DATA which describe\n"
    "               the forward transformation from +orig coordinates\n"
    "               to the 'aaa' coordinate system.\n"
    "             N.B.: The +orig version of this dataset must also be\n"
    "                   readable, since it is also needed when translating\n"
    "                   vectors between SureFit and AFNI coordinates.\n"
    "                   Only the .HEAD files are actually used.\n"
    "\n"
    " -matvec mmm = Read an affine transformation matrix-vector from file\n"
    "               'mmm', which must be in the format\n"
    "                   u11 u12 u13 v1\n"
    "                   u21 u22 u23 v2\n"
    "                   u31 u32 u33 v3\n"
    "               where each 'uij' and 'vi' is a number.  The forward\n"
    "               transformation is defined as\n"
    "                   [ xout ]   [ u11 u12 u13 ] [ xin ]   [ v1 ]\n"
    "                   [ yout ] = [ u21 u22 u23 ] [ yin ] + [ v2 ]\n"
    "                   [ zout ]   [ u31 u32 u33 ] [ zin ]   [ v3 ]\n"
    "\n"
    " Exactly one of -apar or -matvec must be used to specify the\n"
    " transformation.\n"
    "\n"
    " -forward    = -forward means to apply the forward transformation;\n"
    "   *OR*        -backward means to apply the backward transformation\n"
    " -backward     * For example, if the transformation is specified by\n"
    "                  '-apar fred+tlrc', then the forward transformation\n"
    "                  is from +orig to +tlrc coordinates, and the backward\n"
    "                  transformation is from +tlrc to +orig coordinates.\n"
    "               * If the transformation is specified by -matvec, then\n"
    "                  the matrix-vector read in defines the forward\n"
    "                  transform as above, and the backward transformation\n"
    "                  is defined as the inverse.\n"
    "               * If neither -forward nor -backward is given, then\n"
    "                  -forward is the default.\n"
    "\n"
    " -input iii  = Read input 3-vectors from file 'iii' (from stdin if\n"
    "               'iii' is '-' or the -input option is missing).  Input\n"
    "               data may be in one of the following ASCII formats:\n"
    "\n"
    "               * SureFit .coord files:\n"
    "                   BeginHeader\n"
    "                   lines of text ...\n"
    "                   EndHeader\n"
    "                   count\n"
    "                   int x y z\n"
    "                   int x y z\n"
    "                   et cetera...\n"
    "                 In this case, everything up to and including the\n"
    "                 count is simply passed through to the output.  Each\n"
    "                 (x,y,z) triple is transformed, and output with the\n"
    "                 int label that precedes it.  Lines that cannot be\n"
    "                 scanned as 1 int and 3 floats are treated as comments\n"
    "                 and are passed to through to the output unchanged.\n"
    "             N.B.-1: For those using SureFit surfaces created after\n"
    "                     the SureFit/Caret merger (post. 2005), you need\n"
    "                     to use the flag -new_surefit. Talk to Donna about\n"
    "                     this!\n"
    "             N.B.-2: SureFit coordinates are\n"
    "                   x = distance Right    of Left-most      dataset corner\n"
    "                   y = distance Anterior to Posterior-most dataset corner\n"
    "                   z = distance Superior to Inferior-most  dataset corner\n"
    "                 For example, if the transformation is specified by\n"
    "                   -forward -apar fred+tlrc\n"
    "                 then the input (x,y,z) are relative to fred+orig and the\n"
    "                 output (x,y,z) are relative to fred+tlrc.  If instead\n"
    "                   -backward -apar fred+tlrc\n"
    "                 is used, then the input (x,y,z) are relative to fred+tlrc\n"
    "                 and the output (x,y,z) are relative to fred+orig.\n"
    "                 For this to work properly, not only fred+tlrc must be\n"
    "                 readable by Vecwarp, but fred+orig must be as well.\n"
    "                 If the transformation is specified by -matvec, then\n"
    "                 the matrix-vector transformation is applied to the\n"
    "                 (x,y,z) vectors directly, with no coordinate shifting.\n"
    "\n"
    "               * AFNI .1D files with 3 columns\n"
    "                   x y z\n"
    "                   x y z\n"
    "                   et cetera...\n"
    "                 In this case, each (x,y,z) triple is transformed and\n"
    "                 written to the output.  Lines that cannot be scanned\n"
    "                 as 3 floats are treated as comments and are passed\n"
    "                 through to the output unchanged.\n"
    "               N.B.: AFNI (x,y,z) coordinates are in DICOM order:\n"
    "                   -x = Right     +x = Left\n"
    "                   -y = Anterior  +y = Posterior\n"
    "                   -z = Inferior  +z = Superior\n"
    "\n"
    " -output ooo = Write the output to file 'ooo' (to stdout if 'ooo'\n"
    "               is '-', or if the -output option is missing).  If the\n"
    "               file already exists, it will not be overwritten unless\n"
    "               the -force option is also used.\n"
    "\n"
    " -force      = If the output file already exists, -force can be\n"
    "               used to overwrite it.  If you want to use -force,\n"
    "               it must come before -output on the command line.\n"
    "\n"
    "EXAMPLES:\n"
    "\n"
    "  Vecwarp -apar fred+tlrc -input fred.orig.coord > fred.tlrc.coord\n"
    "\n"
    "This transforms the vectors defined in original coordinates to\n"
    "Talairach coordinates, using the transformation previously defined\n"
    "by AFNI markers.\n"
    "\n"
    "  Vecwarp -apar fred+tlrc -input fred.tlrc.coord -backward > fred.test.coord\n"
    "\n"
    "This does the reverse transformation; fred.test.coord should differ from\n"
    "fred.orig.coord only by roundoff error.\n"
    "\n"
    "Author: RWCox - October 2001\n"
   ) ;
   PRINT_COMPILE_DATE ; exit(0) ;
}

/*--------------------------------------------------------------------------*/

static void errex( char * str )
{
   fprintf(stderr,"** FATAL ERROR: %s\n",str) ; exit(1) ;
}

/*--------------------------------------------------------------------------*/

#define NBUF    1024
#define SUREFIT   33
#define AFNI_1D   44

int main( int argc , char *argv[] )
{
   int iarg=1 ;
   FILE *fin=stdin , *fout=stdout ;
   int backward=0 , force=0 ;
   THD_warp *warp=NULL ;
   char lbuf[NBUF] , *cpt ;
   int itype=AFNI_1D , numv=0 , numc=0 ;
   THD_fvec3 vin , vout ;
   float xx,yy,zz ;
   int nn , ii , good=0 ;
   byte old_surefit_mode = 1;
   THD_3dim_dataset *aset=NULL , *oset=NULL ; char *oname=NULL ;

   if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Syntax() ;

   /*-- process command line arguments --*/
   old_surefit_mode = 1;
   while( iarg < argc ){

      /* -input */

      if( strcmp(argv[iarg],"-input") == 0 ){
         if( ++iarg >= argc )
            errex("-input: Need argument after -input") ;
         if( fin != stdin )
            errex("-input: Can't have two -input options") ;
         if( strcmp(argv[iarg],"-") != 0 ){
            fin = fopen( argv[iarg] , "r" ) ;
            if( fin == NULL )
               errex("-input: Can't open input file") ;
         }
         iarg++ ; continue ;
      }

      /* -output */

      if( strcmp(argv[iarg],"-output") == 0 ){
         if( ++iarg >= argc )
            errex("-output: Need argument after -output") ;
         if( fout != stdout )
            errex("-output: Can't have two -output options") ;
         if( strcmp(argv[iarg],"-") != 0 ){
            if( !THD_filename_ok(argv[iarg]) )
               errex("-output: Output filename is illegal") ;
            if( !force && THD_is_file(argv[iarg]) )
               errex("-output: Output file already exists") ;
            fout = fopen( argv[iarg] , "w" ) ;
            if( fout == NULL )
               errex("-output: Can't open output file") ;
         }
         iarg++ ; continue ;
      }

      /* -force */

      if( strcmp(argv[iarg],"-force") == 0 ){
         force = 1 ;
         iarg++ ; continue ;
      }

      /* -new_surefit */

      if( strcmp(argv[iarg],"-new_surefit") == 0 ){
         old_surefit_mode = 0 ;
         iarg++ ; continue ;
      }      
      
      /* -forward */

      if( strcmp(argv[iarg],"-forward") == 0 ){
         backward = 0 ;
         iarg++ ; continue ;
      }

      /* -backward */

      if( strcmp(argv[iarg],"-backward") == 0 ){
         backward = 1 ;
         iarg++ ; continue ;
      }

      /* -apar */

      if( strcmp(argv[iarg],"-apar") == 0 ){
         if( ++iarg >= argc )
            errex("-apar: Need argument after -apar") ;
         if( warp != NULL )
            errex("-apar: Can't specify transformation twice") ;

         /* open dataset with warp */

         aset = THD_open_dataset( argv[iarg] ) ;
         if( !ISVALID_DSET(aset) ){
            sprintf(lbuf,"-apar: Can't open dataset %s\n",argv[iarg]) ;
            errex(lbuf) ;
         }
         if( aset->warp == NULL ){
            sprintf(lbuf,"-apar: Dataset %s does not contain warp",argv[iarg]) ;
            errex(lbuf) ;
         }
         if( aset->view_type == VIEW_ORIGINAL_TYPE ){
            sprintf(lbuf,"-apar: Dataset %s is in the +orig view",argv[iarg]) ;
            errex(lbuf) ;
         }

         /* open +orig version of this dataset */

         sprintf(lbuf,"%s%s+orig.HEAD", aset->dblk->diskptr->directory_name,
                                        aset->dblk->diskptr->prefix         );

         oname = strdup(lbuf) ;
         oset  = THD_open_dataset(lbuf) ;
#if 0
         if( !ISVALID_DSET(oset) ){
            char str[NBUF] ;
            sprintf(str,"-apar: Can't open dataset %s",lbuf) ;
            errex(str) ;
         }
#endif

         warp = aset->warp ;
         iarg++ ; continue ;
      }

      /* -matvec */

      if( strcmp(argv[iarg],"-matvec") == 0 ){
         THD_dvecmat dvm ;
         THD_linear_mapping lmap ;
         
         ZZME(lmap) ;
         
         if( ++iarg >= argc )
            errex("-matvec: Need argument after -matvec") ;
         if( warp != NULL )
            errex("-matvec: Can't specify transformation twice") ;
         dvm = THD_read_dvecmat( argv[iarg] , 0 ) ;
         if( SIZE_DMAT(dvm.mm) == 0.0 )
            errex("-matvec: Can't read transformation") ;

         /* manufacture an AFNI linear map */
         
         lmap.type = MAPPING_LINEAR_TYPE ;
         lmap.mfor.mat[0][0] = dvm.mm.mat[0][0] ;  /* copy matrix in */
         lmap.mfor.mat[0][1] = dvm.mm.mat[0][1] ;
         lmap.mfor.mat[0][2] = dvm.mm.mat[0][2] ;
         lmap.mfor.mat[1][0] = dvm.mm.mat[1][0] ;
         lmap.mfor.mat[1][1] = dvm.mm.mat[1][1] ;
         lmap.mfor.mat[1][2] = dvm.mm.mat[1][2] ;
         lmap.mfor.mat[2][0] = dvm.mm.mat[2][0] ;
         lmap.mfor.mat[2][1] = dvm.mm.mat[2][1] ;
         lmap.mfor.mat[2][2] = dvm.mm.mat[2][2] ;
         lmap.bvec.xyz[0]    = -dvm.vv.xyz[0] ;    /* copy vector in */
         lmap.bvec.xyz[1]    = -dvm.vv.xyz[1] ;
         lmap.bvec.xyz[2]    = -dvm.vv.xyz[2] ;
         LOAD_INVERSE_LMAP(lmap) ;                 /* make inverse map */

         /* manufacture an AFNI warp */

         warp = (THD_warp *) calloc( sizeof(THD_warp) , 1 ) ;
         warp->type = WARP_AFFINE_TYPE ;
         warp->rig_bod.warp = lmap ;

         /* corrected the output location for vv.xyz  17 Aug 2006 [rickr] */
         fprintf(stderr,"++ Using matrix-vector transformation below:\n"
                        "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n"
                        "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n"
                        "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n" ,
       dvm.mm.mat[0][0] , dvm.mm.mat[0][1] , dvm.mm.mat[0][2] , dvm.vv.xyz[0] ,
       dvm.mm.mat[1][0] , dvm.mm.mat[1][1] , dvm.mm.mat[1][2] , dvm.vv.xyz[1] ,
       dvm.mm.mat[2][0] , dvm.mm.mat[2][1] , dvm.mm.mat[2][2] , dvm.vv.xyz[2]) ;

         iarg++ ; continue ;
      }

      /* ERROR */

      { char *str = AFMALL(char, strlen(argv[iarg])+32) ;
        sprintf(str,"Unknown option: %s",argv[iarg]) ;
        errex(str) ;
      }

   } /* end of loop over command line args */

   /*-- check if we are ready --*/

   if( warp == NULL )
      errex("Transformation not specified") ;

   /*-- Read 1st line of input file to determine what to do with it --*/

   cpt = afni_fgets( lbuf , NBUF , fin ) ;
   if( cpt == NULL )
      errex("Couldn't read 1st line from input file") ;

   /*-- check if this is a SureFit file --*/

   if( strstr(cpt,"BeginHeader") != NULL ) itype = SUREFIT ;
   else                                    itype = AFNI_1D ;

   if( itype == SUREFIT && oset == NULL ){
     ERROR_exit(
        "SureFit transformation requires that the +orig dataset also exists:\n"
        "   %s" ,
      oname ) ;
   }

   /*-- if SureFit, echo all header stuff to the output --*/

   if( itype == SUREFIT ){
      int numh=1 ;
      fprintf(fout,"%s",lbuf) ;    /* echo 1st line */
      while(1){                    /* read lines until EndHeader is found */

         cpt = afni_fgets( lbuf , NBUF , fin ) ;
         if( cpt == NULL )
            errex("Input ended before EndHeader was found") ;

         fprintf(fout,"%s",lbuf) ; /* echo line to output */
         numh++ ;

         if( strstr(lbuf,"EndHeader") != NULL ){  /* do next line, too */
            cpt = afni_fgets( lbuf , NBUF , fin ) ;
            if( cpt == NULL )
               errex("Input ended just after EndHeader") ;
            fprintf(fout,"%s",lbuf) ;
            fprintf(stderr,"++ Wrote %d SureFit header lines to output\n",numh);
            break ;                /* end of loop */
         }
      }

      cpt = afni_fgets( lbuf , NBUF , fin ) ;  /* get next line, with 1st vector */
      if( cpt == NULL )
         errex("Input ended just after Node count") ;

   } /* end of SureFit header echo */

   /*-- From this point on, process each line with 1 vector. --*/
   /*-- At the start of the loop, 1 line's data is in lbuf.  --*/

   do{

      switch( itype ){
         case SUREFIT:
            ii = sscanf(lbuf,"%d%f%f%f",&nn,&xx,&yy,&zz) ;
            good = (ii == 4) ;
         break ;

         case AFNI_1D:
            ii = sscanf(lbuf,"%f%f%f",&xx,&yy,&zz) ;
            good = (ii == 3) ;
         break ;
      }

      if( !good ){  /* just echo line to output */

         fprintf(fout,"%s",lbuf) ;
         numc++ ;

      } else {      /* transform a vector!!! */

         LOAD_FVEC3(vin,xx,yy,zz) ;
#ifdef DEBUG
fprintf(stderr,"\n") ;
DUMP_FVEC3("vin              ",vin) ;
#endif

         if( itype == SUREFIT && old_surefit_mode && aset != NULL ){
            if( backward ) vin = THD_surefit_to_dicomm( aset , vin ) ;
            else           vin = THD_surefit_to_dicomm( oset , vin ) ;
#ifdef DEBUG
DUMP_FVEC3("surefit_to_dicomm",vin) ;
#endif
         }

         if( backward ) vout = myAFNI_backward_warp_vector( warp , vin ) ;
         else           vout = myAFNI_forward_warp_vector ( warp , vin ) ;

#ifdef DEBUG
DUMP_FVEC3("vout             ",vout) ;
#endif

         if( itype == SUREFIT && old_surefit_mode && aset != NULL ){
            if( backward ) vout = THD_dicomm_to_surefit( oset , vout ) ;
            else           vout = THD_dicomm_to_surefit( aset , vout ) ;
#ifdef DEBUG
DUMP_FVEC3("dicomm_to_surefit",vout) ;
#endif
         }

         if( itype == SUREFIT ) fprintf(fout,"%d ",nn) ;
         fprintf(fout,"%f %f %f\n",vout.xyz[0],vout.xyz[1],vout.xyz[2]) ;
         numv++ ;
      }

      cpt = afni_fgets( lbuf , NBUF , fin ) ;  /* get next line */

   } while( cpt != NULL ) ;  /* loop until no data can be read */

   if( numc > 0 )
      fprintf(stderr,"++ Wrote %d vectors;  %d comments\n",numv,numc) ;
   else
      fprintf(stderr,"++ Wrote %d vectors\n",numv) ;

   exit(0) ;
}

/*--------------------------------------------------------------------------
  The following functions are adapted from afni.c
----------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
   Forward transform a vector following a warp
--------------------------------------------------------------------------*/

THD_fvec3 myAFNI_forward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
{
   THD_fvec3 new_fv ;

   if( warp == NULL ) return old_fv ;

   switch( warp->type ){

      default: new_fv = old_fv ; break ;

      case WARP_TALAIRACH_12_TYPE:{
         THD_linear_mapping map ;
         int iw ;

         /* forward transform each possible case,
            and test if result is in bot..top of defined map */

         for( iw=0 ; iw < 12 ; iw++ ){
            map    = warp->tal_12.warp[iw] ;
            new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;

            if( new_fv.xyz[0] >= map.bot.xyz[0] &&
                new_fv.xyz[1] >= map.bot.xyz[1] &&
                new_fv.xyz[2] >= map.bot.xyz[2] &&
                new_fv.xyz[0] <= map.top.xyz[0] &&
                new_fv.xyz[1] <= map.top.xyz[1] &&
                new_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
         }
      }
      break ;

      case WARP_AFFINE_TYPE:{
         THD_linear_mapping map = warp->rig_bod.warp ;
         new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
      }
      break ;

   }
   return new_fv ;
}

/*------------------------------------------------------------------------
   Backward transform a vector following a warp
--------------------------------------------------------------------------*/

THD_fvec3 myAFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
{
   THD_fvec3 new_fv ;

   if( warp == NULL ) return old_fv ;

   switch( warp->type ){

      default: new_fv = old_fv ; break ;

      case WARP_TALAIRACH_12_TYPE:{
         THD_linear_mapping map ;
         int iw ;

         /* test if input is in bot..top of each defined map */

         for( iw=0 ; iw < 12 ; iw++ ){
            map = warp->tal_12.warp[iw] ;

            if( old_fv.xyz[0] >= map.bot.xyz[0] &&
                old_fv.xyz[1] >= map.bot.xyz[1] &&
                old_fv.xyz[2] >= map.bot.xyz[2] &&
                old_fv.xyz[0] <= map.top.xyz[0] &&
                old_fv.xyz[1] <= map.top.xyz[1] &&
                old_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
         }
         new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
      }
      break ;

      case WARP_AFFINE_TYPE:{
         THD_linear_mapping map = warp->rig_bod.warp ;
         new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
      }
      break ;

   }
   return new_fv ;
}

#if 0
/*--------------------------------------------------------------------------
   The following routines are used to convert DICOM order coordinates
   (used in AFNI) to SureFit order coordinates
----------------------------------------------------------------------------*/

static THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *dset , THD_fvec3 fv )
{
   float xx,yy,zz , xbase,ybase,zbase ;
   THD_fvec3 vout ;
#ifdef DEBUG
static int first=1 ;
#endif

   xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;   /* xyz now LPI */

   if( dset != NULL ){
      THD_fvec3 v1 , v2 ;
      LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
      v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
      LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
                     DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
                     DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset)  ) ;
      v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
      xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;  /* Left-most */
      ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;  /* Posterior */
      zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;                   /* Inferior  */
   } else {
      xbase = ybase = zbase = 0.0 ;
   }
#ifdef DEBUG
if(first){fprintf(stderr,"d2s base=%f %f %f\n",xbase,ybase,zbase);first=0;}
#endif

   vout.xyz[0] = xx - xbase ;
   vout.xyz[1] = yy - ybase ;
   vout.xyz[2] = zz - zbase ; return vout ;
}

static THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 fv )
{
   float xx,yy,zz , xbase,ybase,zbase ;
   THD_fvec3 vout ;
#ifdef DEBUG
static int first=1 ;
#endif

   xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;   /* xyz now RAI */

   if( dset != NULL ){
      THD_fvec3 v1 , v2 ;
      LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
      v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
      LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
                     DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
                     DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset)  ) ;
      v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
      xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
      ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
      zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
   } else {
      xbase = ybase = zbase = 0.0 ;
   }
#ifdef DEBUG
if(first){fprintf(stderr,"s2d base=%f %f %f\n",xbase,ybase,zbase);first=0;}
#endif

   vout.xyz[0] = xx - xbase ;
   vout.xyz[1] = yy - ybase ;
   vout.xyz[2] = zz + zbase ; return vout ;
}
#endif
