/*------------------------------------------------------------------------
     ***  This program does something, but nobody is sure what.  ***
     ***  14 Apr 2003 - RWCox.                                   ***
--------------------------------------------------------------------------*/

/*

  [PT: Sep 9, 2019] add in opt: -disp_obl_xform_only
                    + output transform between oblique coords.

*/

#include "mrilib.h"

/*--------------------------------------------------------------------------*/
#undef DEBUG_ON

#define MATVEC_FOR 1
#define MATVEC_BAC 2

static THD_vecmat dicom_in2out , dicom_out2in ;  /* coordinate warps */
/*float compute_oblique_angle(mat44 ijk_to_dicom44);*/
static void Compute_Deoblique_Transformation(THD_3dim_dataset *dset,mat44 *Tw);

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

void warp_dicom_in2out( float  xin , float  yin , float  zin ,
                        float *xout, float *yout, float *zout )
{
   THD_fvec3 xxx ;
   LOAD_FVEC3( xxx , xin,yin,zin ) ;
   xxx = VECMAT_VEC( dicom_in2out , xxx ) ;
   UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
}

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

void warp_dicom_out2in( float  xin , float  yin , float  zin ,
                        float *xout, float *yout, float *zout )
{
   THD_fvec3 xxx ;
   LOAD_FVEC3( xxx , xin,yin,zin ) ;
   xxx = VECMAT_VEC( dicom_out2in , xxx ) ;
   UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
}

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

int main( int argc , char * argv[] )
{
   THD_3dim_dataset *inset , *outset=NULL , *newgset=NULL ;
   char * prefix = "warped" ;
   int nopt=1 , verb=0 , zpad=0 , fsl=0 , gsetopt=0;
   int use_matvec=0 ;
   int use_newgrid=0 ;
   float ddd_newgrid=0.0 ;
   void *newggg=NULL ; int gflag=0 ;

   int tta2mni=0 , mni2tta=0 ;   /* 11 Mar 2004 */
   int matdefined=0 ;

   THD_coorder cord ;
   ATR_float *atr_matfor=NULL , *atr_matinv=NULL ;
   THD_dmat33 tmat ;
   THD_dfvec3 tvec ;
   mat44 Tw, Tc, Tw_inv, Tr, Tw2;
   int oblique_flag = 0;
   THD_3dim_dataset *oblparset=NULL ;
   float angle; 
   float matar[12] ; char anam[64] ;
#if 0
   MRI_IMAGE *matflim=NULL ;
   float     *matflar=NULL ;
#endif

   int DISP_OBL_XFORM_ONLY = 0 ;  // [PT: Sep 13, 2019] for getting oblique
                            // transform matrix most easily

   /*-- help? --*/

   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
     printf("\n"
            "Usage: 3dWarp [options] dataset\n"
            "Warp (spatially transform) one 3D dataset.\n"
            "\n"
            "* Here, the 'warp' is a 3x4 matrix = affine transform of space which\n"
            "  the user supplies.\n"
            "* Such a transformation can be computed by 3dWarpDrive or 3dAllineate,\n"
            "  by matching one dataset to another.\n"
            "* However, 3dAllineate can also transform a dataset given a matrix,\n"
            "  so the usefulness of these older 3dWarp programs is now limited.\n"
            "*****\n"
            "***** For nonlinear spatial warping, see program 3dQwarp. *****\n"
            "*****\n"
            "--------------------------\n"
            "Transform Defining Options: [exactly one of these must be used]\n"
            "--------------------------\n"
            "  -matvec_in2out mmm = Read a 3x4 affine transform matrix+vector\n"
            "                        from file 'mmm':\n"
            "                         x_out = Matrix x_in + Vector\n"
            "\n"
            "  -matvec_out2in mmm = Read a 3x4 affine transform matrix+vector\n"
            "                         from file 'mmm':\n"
            "                         x_in = Matrix x_out + Vector\n"
            "\n"
            "     ** N.B.: The coordinate vectors described above are\n"
            "               defined in DICOM ('RAI') coordinate order.\n"
            "               (Also see the '-fsl_matvec option, below.)\n"
            "     ** N.B.: Using the special name 'IDENTITY' for 'mmm'\n"
            "               means to use the identity matrix.\n"
            "     ** N.B.: You can put the matrix on the command line\n"
            "               directly by using an argument of the form\n"
            "       'MATRIX(a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34)'\n"
            "               in place of 'mmm', where the aij values are the\n"
            "               matrix entries (aij = i-th row, j-th column),\n"
            "               separated by commas.\n"
            "             * You will need the 'forward single quotes' around\n"
            "               the argument.\n"
            "\n"
            "  -tta2mni = Transform a dataset in Talairach-Tournoux Atlas\n"
            "              coordinates to MNI-152 coordinates.\n"
            "  -mni2tta = Transform a dataset in MNI-152 coordinates to\n"
            "              Talairach-Tournoux Atlas coordinates.\n"
            "\n"
            "  -matparent mset = Read in the matrix from WARPDRIVE_MATVEC_*\n"
            "                     attributes in the header of dataset 'mset',\n"
            "                     which must have been created by program\n"
            "                     3dWarpDrive.  In this way, you can apply\n"
            "                     a transformation matrix computed from\n"
            "                     in 3dWarpDrive to another dataset.\n\n"
            "     ** N.B.: The above option is analogous to the -rotparent\n"
            "                option in program 3drotate.  Use of -matparent\n"
            "                should be limited to datasets whose spatial\n"
            "                coordinate system corresponds to that which\n"
            "                was used for input to 3dWarpDrive (i.e., the\n"
            "                input to 3dWarp should overlay properly with\n"
            "                the input to 3dWarpDrive that generated the\n"
            "                -matparent dataset).\n\n"
  "  -card2oblique   obl_dset \n"
  "       or \n"
  "  -oblique_parent obl_dset = Read in the oblique transformation matrix\n"
  "     from an oblique dataset and make cardinal dataset oblique to match.\n"
  "  -deoblique or\n"
  "  -oblique2card = Transform an oblique dataset to a cardinal dataset\n"
  "     Both these oblique transformation options require a new grid for the\n"
  "     output as specified with the -newgrid or -gridset options\n"
  "     or a new grid will be assigned based on the minimum voxel spacing\n"
  "    ** N.B.: EPI time series data should be time shifted with 3dTshift before"
  "                rotating the volumes to a cardinal direction\n"
  "\n"
  "  -disp_obl_xform_only = (new opt) just display the obliquity transform\n"
  "                         matrix that would be applied to make the output\n"
  "                         dset;  very useful for moving between oblique\n"
  "                         coords, such as with '-oblique_parent ..' or\n"
  "                         '-deoblique'.  The result can be dumped into a\n"
  "                         text file, e.g.:\n" 
  "                              > textfile.aff12.1D\n"
  "                         No dataset is created or changed.\n"
  "\n"
  "Sample usages:\n"
  " 3dWarpDrive -affine_general -base d1+orig -prefix d2WW -twopass"
    " -input d2+orig\n"
  " 3dWarp -matparent d2WW+orig -prefix epi2WW epi2+orig\n\n"
  " 3dWarp -card2oblique oblique_epi+orig -prefix oblique_anat card_anat+orig\n"
  " 3dWarp -oblique2card -prefix card_epi_tshift -newgrid 3.5 epi_tshift+orig\n"
  " 3dWarp -card2oblique oblique_epi.nii -disp_obl_transf_only epi_tshift+orig\\\n"
  "        > mat_obl_transform.aff12.1D\n"
  "\n"
  "Example of warping +tlrc results back to +orig space of some subject\n"
  "(get xform matrix, apply it, tell dataset it is not in orig space):\n"
  "\n"
  "    cat_matvec subj1_anat+tlrc::WARP_DATA > tlrc_xform.1D\n"
  "    3dWarp -matvec_out2in tlrc_xform.1D -prefix group_warped+tlrc \\\n"
  "           -gridset subj1_epi+orig -cubic group_data+tlrc\n"
  "    3drefit -view orig group_warped+tlrc\n"
#if 0
  "\n"
  "  -matfile mname  = Read in the file 'mname', which consists\n"
  "                     of 12 ASCII-formatted numbers per line.\n"
  "                     The i-th input line defines a 3x4 matrix\n"
  "                     which is used to transform the i-th\n"
  "                     sub-brick of the input dataset.\n"
#endif
  "\n"
  "-----------------------\n"
  "Other Transform Options:\n"
  "-----------------------\n"
  "  -linear     }\n"
  "  -cubic      } = Chooses spatial interpolation method.\n"
  "  -NN         } =   [default = linear]\n"
  "  -quintic    }\n"
  "  -wsinc5     }\n"
  "\n"
  "  -fsl_matvec   = Indicates that the matrix file 'mmm' uses FSL\n"
  "                    ordered coordinates ('LPI').  For use with\n"
  "                    matrix files from FSL and SPM.\n"
  "\n"
  "  -newgrid ddd  = Tells program to compute new dataset on a\n"
  "                    new 3D grid, with spacing of 'ddd' mmm.\n"
  "                  * If this option is given, then the new\n"
  "                    3D region of space covered by the grid\n"
  "                    is computed by warping the 8 corners of\n"
  "                    the input dataset, then laying down a\n"
  "                    regular grid with spacing 'ddd'.\n"
  "                  * If this option is NOT given, then the\n"
  "                    new dataset is computed on the old\n"
  "                    dataset's grid.\n"
  "\n"
  "  -gridset ggg  = Tells program to compute new dataset on the\n"
  "                    same grid as dataset 'ggg'.\n"
  "\n"
  "  -zpad N       = Tells program to pad input dataset with 'N'\n"
  "                    planes of zeros on all sides before doing\n"
  "                    transformation.\n"
  "---------------------\n"
  "Miscellaneous Options:\n"
  "---------------------\n"
  "  -verb         = Print out some information along the way.\n"
  "  -prefix ppp   = Sets the prefix of the output dataset.\n"
  "\n"
  ) ;
     PRINT_COMPILE_DATE ; exit(0) ;
   }

   /*-- startup mechanics --*/

   mainENTRY("3dWarp main"); machdep(); AFNI_logger("3dWarp",argc,argv);
   PRINT_VERSION("3dWarp") ; AUTHOR("RW Cox") ;

   THD_coorder_fill( my_getenv("AFNI_ORIENT") , &cord ) ;  /* 12 Mar 2004 */

   /*-- command line options --*/

   while( nopt < argc && argv[nopt][0] == '-' ){

     /*----- 06 Jul 2005 -----*/

#if 0
     if( strncmp(argv[nopt],"-matfile",8) == 0 ){
       if( matdefined )
         ERROR_exit("-matfile: Matrix already defined!\n");
       if( ++nopt >= argc )
         ERROR_exit("Need argument after -matfile!\n");
       matflim = mri_read_ascii( argv[nopt] ) ;
       if( matflim == NULL )
         ERROR_exit("Can't read from -matfile %s\n",argv[nopt]);
       if( matflim->nx < 12 )
         ERROR_exit("Need 12 numbers per line in -matfile %s\n",argv[nopt]) ;

       matflar = MRI_FLOAT_PTR(matflim) ;
       matdefined = 1 ; nopt++ ; continue ;
     }
#endif

     /*-----*/

     if( strncmp(argv[nopt],"-matparent",10) == 0 ){
       ATR_float *atr ; float *matar ; THD_3dim_dataset *matparset=NULL ;

       if( matdefined )
         ERROR_exit("-matparent: Matrix already defined!\n");
       if( ++nopt >= argc )
         ERROR_exit("Need argument after -matparent!\n");
       matparset = THD_open_dataset( argv[nopt] ) ;
       if( matparset == NULL )
         ERROR_exit("Can't open -matparent %s\n",argv[nopt]);
       atr = THD_find_float_atr( matparset->dblk, 
                                 "WARPDRIVE_MATVEC_INV_000000" );
       if( atr != NULL ) 
         atr_matinv = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
       atr = THD_find_float_atr( matparset->dblk, 
                                 "WARPDRIVE_MATVEC_FOR_000000" );
       if( atr != NULL ) 
         atr_matfor = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
       if( atr_matinv == NULL ||  atr_matinv->nfl < 12 )
         ERROR_exit( "-matparent %s doesn't have WARPDRIVE attributes!?\n",
                     argv[nopt]) ;

       matar = atr_matinv->fl ;
       if( strstr(argv[nopt-1],"INV") == NULL ){
         LOAD_MAT  ( dicom_in2out.mm, matar[0],matar[1],matar[2],
                                      matar[4],matar[5],matar[6],
                                      matar[8],matar[9],matar[10] ) ;
         LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;
         dicom_out2in = INV_VECMAT( dicom_in2out ) ;
       } else {
         LOAD_MAT  ( dicom_out2in.mm, matar[0],matar[1],matar[2],
                                      matar[4],matar[5],matar[6],
                                      matar[8],matar[9],matar[10] ) ;
         LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;
         dicom_in2out = INV_VECMAT( dicom_out2in ) ;
       }

       DSET_delete(matparset) ;

       use_matvec = MATVEC_FOR ; matdefined = 1 ; nopt++ ; continue ;
     }

     /*----- 11 Mar 2004 -----*/

     if( strcmp(argv[nopt],"-tta2mni") == 0 ){
       if( matdefined )
         ERROR_exit("-tta2mni: Matrix already defined!\n");
       matdefined = 1 ; tta2mni = 1 ; nopt++ ; continue ;
     }

     if( strcmp(argv[nopt],"-mni2tta") == 0 ){
       if( matdefined )
         ERROR_exit("-mni2tta: Matrix already defined!\n");
       matdefined = 1 ; mni2tta = 1 ; nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-fsl_matvec") == 0 ){
       fsl = 1 ; nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-gridset") == 0 ){
       if( use_newgrid )
         ERROR_exit("Can't use -gridset with -newgrid!\n");
       if( newgset != NULL )
         ERROR_exit("Can't use -gridset twice!\n");
       if( ++nopt >= argc )
         ERROR_exit("Need argument after -gridset!\n");
       newgset = (THD_3dim_dataset *) -1; /* open the dataset later instead */
       gsetopt = nopt;
       nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-zpad") == 0 ){
       if( ++nopt >= argc )
         ERROR_exit("Need argument after -zpad!\n");
       zpad = (int) strtod( argv[nopt] , NULL ) ;
       if( zpad < 0 )
         ERROR_exit("Illegal negative argument after -zpad!\n");
       else if( zpad == 0 )
         INFO_message("NOTA BENE: -zpad is 0 ==> ignored\n") ;
       nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-newgrid") == 0 ){
       if( newgset != NULL )
         ERROR_exit("Can't use -newgrid with -gridset!\n");
       if( use_newgrid )
         ERROR_exit("Can't use -newgrid twice!\n");
       if( ++nopt >= argc )
         ERROR_exit("Need argument after -newgrid!\n");

       ddd_newgrid = strtod( argv[nopt] , NULL ) ;
       if( ddd_newgrid <= 0.0 )
         ERROR_exit("Illegal argument after -newgrid!\n");

       use_newgrid = 1 ; nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-NN")     == 0 ){
       mri_warp3D_method( MRI_NN )     ; nopt++ ; continue ;
     }
     if( strcmp(argv[nopt],"-linear") == 0 ){
       mri_warp3D_method( MRI_LINEAR ) ; nopt++ ; continue ;
     }
     if( strcmp(argv[nopt],"-cubic")  == 0 ){
       mri_warp3D_method( MRI_CUBIC )  ; nopt++ ; continue ;
     }
     if( strcmp(argv[nopt],"-quintic") == 0 ){
       mri_warp3D_method( MRI_QUINTIC ); nopt++ ; continue ;  /* 06 Aug 2003 */
     }

     /* rcr - todo: move e.g. INTERP_methname to libmri and apply as generic */
     if( strcmp(argv[nopt],"-wsinc5") == 0 ){  /* 8 Feb 2019 [rickr] */
       mri_warp3D_method( MRI_WSINC5 ); nopt++ ; continue ;
     }

     /*-----*/

     if( strcmp(argv[nopt],"-prefix") == 0 ){
       if( ++nopt >= argc )
         ERROR_exit("Need an argument after -prefix!\n");
       if( !THD_filename_ok(argv[nopt]) )
         ERROR_exit("-prefix argument is invalid!\n");

       prefix = argv[nopt] ; nopt++ ; continue ;
     }

     /*-----*/

     if( strncmp(argv[nopt],"-verbose",5) == 0 ){
       verb++ ; nopt++ ; continue ;
     }

     /*-----*/

     if( strncmp(argv[nopt],"-matvec_",8) == 0 ){
       MRI_IMAGE *matim ; float *matar , dm ;

       if( matdefined )
         ERROR_exit("%s: matrix already defined!\n",argv[nopt]);
       if( ++nopt >= argc )
         ERROR_exit("Need an argument after -matvec!\n");
       if( strcmp(argv[nopt],"IDENTITY") == 0 ){ /* load identity matrix */
         matim = mri_new( 4,3 , MRI_float ) ;    /* will be all zero */
         matar = MRI_FLOAT_PTR(matim) ;
         matar[0] = matar[5] = matar[10] = 1.0 ;
       } else if( strncmp(argv[nopt],"MATRIX(",7) == 0 ){ /* matrix from arg */
         int nn ;
         matim = mri_new( 4,3 , MRI_float ) ;    /* will be all zero */
         matar = MRI_FLOAT_PTR(matim) ;
         nn = sscanf(argv[nopt],"MATRIX(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f",
                     matar+0 , matar+1 , matar+2 , matar+3 ,
                     matar+4 , matar+5 , matar+6 , matar+7 ,
                     matar+8 , matar+9 , matar+10, matar+11 ) ;
         if( nn < 12 )
           ERROR_exit("Could only decode %d numbers from %s\n",nn,argv[nopt]) ;
       } else {
         matim = mri_read_ascii( argv[nopt] ) ;
         if( matim == NULL )
           ERROR_exit("Can't read -matvec file %s\n",argv[nopt]);
         if( matim->nx != 4 || matim->ny < 3 )
           ERROR_exit("-matvec file not a 3x4 matrix!\n");
         matar = MRI_FLOAT_PTR(matim) ;
       }
       
       if (!strcmp(argv[nopt-1],"-matvec_in2out")) {
         use_matvec = MATVEC_FOR;  
       } else if (!strcmp(argv[nopt-1],"-matvec_out2in")) {
         use_matvec = MATVEC_BAC;  
       } else {
         ERROR_exit( "Your -matvec option must end with either \n"
                     "_in2out or _out2in. You have %s\n",
                     argv[nopt-1]+7);
       }
       switch( use_matvec ){

         case MATVEC_FOR:
           LOAD_MAT  ( dicom_in2out.mm, matar[0],matar[1],matar[2],
                                        matar[4],matar[5],matar[6],
                                        matar[8],matar[9],matar[10] ) ;
           LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;

           dm = MAT_DET( dicom_in2out.mm ) ;
           if( dm == 0.0 )
             ERROR_exit("Determinant of matrix is 0\n");

           dicom_out2in = INV_VECMAT( dicom_in2out ) ;
         break ;

         case MATVEC_BAC:
           LOAD_MAT  ( dicom_out2in.mm, matar[0],matar[1],matar[2],
                                        matar[4],matar[5],matar[6],
                                        matar[8],matar[9],matar[10] ) ;
           LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;

           dm = MAT_DET( dicom_out2in.mm ) ;
           if( dm == 0.0 )
             ERROR_exit("Determinant of matrix is 0\n");

           dicom_in2out = INV_VECMAT( dicom_out2in ) ;
         break ;
       }

       matdefined = 1 ; nopt++ ; continue ;
     }

     if((strncmp(argv[nopt],"-oblique_parent",15) == 0 ) || \
        (strncmp(argv[nopt],"-card2oblique",13) == 0 )) {
       if( matdefined )
         ERROR_exit("-oblique_parent: Matrix already defined!\n");
      THD_set_oblique_report(0,0); /* turn off obliquity warning */
      if( ++nopt >= argc )
         ERROR_exit("Need argument after -oblique_parent!\n");
       oblparset = THD_open_dataset( argv[nopt] ) ;
       if( oblparset == NULL )
         ERROR_exit("Can't open -oblique_parent %s\n",argv[nopt]);
       oblique_flag = 2;  matdefined = 1; nopt++ ; continue ;
     }

     if((strncmp(argv[nopt],"-deoblique",10) == 0 ) ||
       ( strncmp(argv[nopt],"-oblique2card",13) == 0 )){
        THD_set_oblique_report(0,0);
        oblique_flag = 1; matdefined = 1; nopt++ ; continue ;
     }

     // [PT: Sep 13, 2019] show transform and exit
     if((strncmp(argv[nopt],"-disp_obl_xform_only",20) == 0 )){
        DISP_OBL_XFORM_ONLY = 1; nopt++ ; continue ;
     }

     /*-----*/

     ERROR_message("Unknown option %s\n",argv[nopt]) ;
     suggest_best_prog_option(argv[0], argv[nopt]);
     exit(1);
   } /* end of loop over options */

   // ===============================================================

   /*-- check to see if we have a warp specified --*/

   if( !matdefined )
     ERROR_exit("No transformation on command line!?\n");

   if( tta2mni || mni2tta ) fsl = 0 ;

   /*-- 06 Aug 2003: if FSL matrix, flip stuff around --*/

   if( use_matvec && fsl ){
     dicom_in2out.mm.mat[0][2] = -dicom_in2out.mm.mat[0][2] ;
     dicom_in2out.mm.mat[1][2] = -dicom_in2out.mm.mat[1][2] ;
     dicom_in2out.mm.mat[2][0] = -dicom_in2out.mm.mat[2][0] ;
     dicom_in2out.mm.mat[2][1] = -dicom_in2out.mm.mat[2][1] ;
     dicom_in2out.vv.xyz[0]    = -dicom_in2out.vv.xyz[0]    ;
     dicom_in2out.vv.xyz[1]    = -dicom_in2out.vv.xyz[1]    ;

     dicom_out2in.mm.mat[0][2] = -dicom_out2in.mm.mat[0][2] ;
     dicom_out2in.mm.mat[1][2] = -dicom_out2in.mm.mat[1][2] ;
     dicom_out2in.mm.mat[2][0] = -dicom_out2in.mm.mat[2][0] ;
     dicom_out2in.mm.mat[2][1] = -dicom_out2in.mm.mat[2][1] ;
     dicom_out2in.vv.xyz[0]    = -dicom_out2in.vv.xyz[0]    ;
     dicom_out2in.vv.xyz[1]    = -dicom_out2in.vv.xyz[1]    ;
   }

   /*-- 1 remaining argument should be a dataset --*/

   if( nopt != argc-1 ){
     ERROR_message("Command line should have exactly 1 dataset!") ;
     ERROR_exit   ("Whereas there seems to be %d of them!\n", argc-nopt ) ;
   }

   /*-- input dataset header --*/

   inset = THD_open_dataset( argv[nopt] ) ;
   if( !ISVALID_DSET(inset) )
     ERROR_exit("Can't open dataset %s\n",argv[nopt]);

   #if 0 
   {
      /*- 25 Mar 2008: Deal with xforms that apply to i_in rather 
                       than x_in 
          WARNING: This code has not been tested for lack of data.
          Say the xform specified on command line is M_i such that
               x_out = M_i i_in +  Vector,
               where i_in is a voxel's IJK index, rather than a voxel's
               dicomm coordinate. 
               or 
               X_out = M_I I_in (using 4x4 xform matrix convention) 
          
          We want to calculate the transform M_X such that:
               X_out = M_X X_in , 
               where X_in is the voxel's dicomm coordinate. 
          
               So M_X X_in = M_I I_in
          
          In AFNI's dset structure,
               X_in = M_ijk I_in , 
               where M_ijk is daxes->ijk_to_dicom_real
          
          Now we have:
               M_X M_ijk I_in = M_I I_in
               or
               M_X = M_I inv(M_ijk)
          
          So if needed, the block below will replace
          dicom_in2out with M_X . But first we have
          to wait for some data.         
          */
      fprintf(stderr,"WARNING\nWARNING!\nWARNING!!\n");
      if(ISVALID_MAT44(inset->daxes->ijk_to_dicom_real)) {
         /* load dicom_in2out into Tw */
         LOAD_MAT44(Tw, 
          dicom_in2out.mm.mat[0][0], 
          dicom_in2out.mm.mat[0][1], 
          dicom_in2out.mm.mat[0][2], dicom_in2out.vv.xyz[0],
          dicom_in2out.mm.mat[1][0], 
          dicom_in2out.mm.mat[1][1], 
          dicom_in2out.mm.mat[1][2], dicom_in2out.vv.xyz[1],
          dicom_in2out.mm.mat[2][0], 
          dicom_in2out.mm.mat[2][1], 
          dicom_in2out.mm.mat[2][2], dicom_in2out.vv.xyz[2]);
         DUMP_MAT44("MI:\n", Tw);
         /* equivalent dicom transform is Tw*inv(ijk_to_dicom_real) */
         Tc = MAT44_INV(inset->daxes->ijk_to_dicom_real);
         DUMP_MAT44("inv(Mijk):\n", Tc);
         Tw2 = MAT44_MUL(Tw,Tc);
         DUMP_MAT44("MI*inv(Mijk):\n", Tw2);
         /* Now reload Tw2 into dicom_in2out */
         LOAD_MAT  ( dicom_in2out.mm, 
                     Tw2.m[0][0],Tw2.m[0][1],Tw2.m[0][2],
                     Tw2.m[1][0],Tw2.m[1][1],Tw2.m[1][2],
                     Tw2.m[2][0],Tw2.m[2][1],Tw2.m[2][2] ) ;
         LOAD_FVEC3( dicom_in2out.vv, 
                     Tw2.m[0][3],Tw2.m[1][3],Tw2.m[2][3] ) ;
         /* And recalculate the inverse of dicom_in2out */
         dicom_out2in = INV_VECMAT( dicom_in2out ) ;
      }
   }
   #endif
/* these tests don't fit in well with other changes for template spaces */
#if 0
   /*- 11 Mar 2004: check for coherency -*/

   if( tta2mni || mni2tta ){
     if( inset->view_type != VIEW_TALAIRACH_TYPE ){
       WARNING_message("Input dataset not in +tlrc view!\n") ;
     } else {
       ATR_string *ast ;
       ast = THD_find_string_atr( inset->dblk , "TLRC_SPACE" ) ;
       if( ast != NULL ){
         if( strcmp(ast->ch,"MNI-152") == 0 && tta2mni ){
           WARNING_message("Input dataset appears to be MNI-152 already!\n");
         }
       }
     }
   }
#endif

   /*-- do the heavy lifting --*/

   if( use_newgrid ){
     newggg = &ddd_newgrid ; gflag = WARP3D_NEWGRID ;
   } else if( newgset != NULL ){
     newggg = newgset = THD_open_dataset( argv[gsetopt] ) ;
       if( newgset == NULL )
         ERROR_exit("Can't open -gridset %s\n",argv[nopt]);
     gflag = WARP3D_NEWDSET ;
   }

    /* handling oblique and deoblique cases */
   if(oblique_flag) {
      float *matar, dm;
      if(oblique_flag==1)   /* if deobliquing, everything is in the same dataset */
         oblparset = inset;
      Compute_Deoblique_Transformation(oblparset, &Tw);

      if(oblique_flag==1)   /* deoblique case*/
         matar = &Tw.m[0][0];
      else {               /* obliquing case */
         /* don't need oblique parent dataset anymore */
         DSET_delete(oblparset) ;
         Tw_inv = MAT44_INV(Tw);
         matar = &Tw_inv.m[0][0];
         Tw = Tw_inv;      /* just for the DUMP */
#ifdef DEBUG_ON
DUMP_MAT44("Tw_inv",Tw_inv);
#endif
	 if(ISVALID_MAT44(inset->daxes->ijk_to_dicom_real)) {
	   angle = THD_compute_oblique_angle(inset->daxes->ijk_to_dicom_real, 0);
	   if(angle>0.0) {
              INFO_message("Deobliquing original dataset before obliquing\n");
              INFO_message("  Combining oblique transformations");
              Compute_Deoblique_Transformation(inset, &Tw2);
              Tw = MAT44_MUL(Tw_inv, Tw2);
#ifdef DEBUG_ON
DUMP_MAT44("Twcombined", Tw);
#endif
              matar = &Tw.m[0][0];
           }
         }
      }

      /* show overall oblique/deoblique transformation */
      if(oblique_flag)
         if(DISP_OBL_XFORM_ONLY || oblique_flag )
            DUMP_MAT44("Obliquity Transformation :", Tw);
      
      if( DISP_OBL_XFORM_ONLY )
         exit(0);

      LOAD_MAT  ( dicom_in2out.mm, matar[0],matar[1],matar[2],
                                   matar[4],matar[5],matar[6],
                                   matar[8],matar[9],matar[10] ) ;
      LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;

      dm = MAT_DET( dicom_in2out.mm ) ;
      if( dm == 0.0 )
        ERROR_exit("Determinant of matrix is 0\n");

      dicom_out2in = INV_VECMAT( dicom_in2out ) ;

      use_matvec = MATVEC_FOR ; matdefined = 1 ;
      if(newggg==NULL) {
         ddd_newgrid = DSET_MIN_DEL(inset);
         newggg = &ddd_newgrid ; gflag = WARP3D_NEWGRID ;
         INFO_message("Using minimum spacing of %f mm for new grid spacing",ddd_newgrid);
      }
    }
   /* if data is not being deobliqued or obliquified */
   if(!oblique_flag) {
     if(ISVALID_MAT44(inset->daxes->ijk_to_dicom_real)) {
       angle = THD_compute_oblique_angle(inset->daxes->ijk_to_dicom_real,1);
       if(angle>0.0) {  
         THD_dicom_card_xform(inset, &tmat, &tvec); 
         LOAD_MAT44(Tc, 
          tmat.mat[0][0], tmat.mat[0][1], tmat.mat[0][2], tvec.xyz[0],
          tmat.mat[1][0], tmat.mat[1][1], tmat.mat[1][2], tvec.xyz[1],
          tmat.mat[2][0], tmat.mat[2][1], tmat.mat[2][2], tvec.xyz[2]);
         Tr = MAT44_SUB(Tc,inset->daxes->ijk_to_dicom_real);
         if(MAT44_NORM(Tr)>0.001) {
            WARNING_message("Deoblique datasets with 3dWarp before proceeding"
               " with other transformations");
         }
       }
      }
   }

   if( use_matvec ){
     outset = THD_warp3D_affine( inset , dicom_out2in ,
                                 newggg , prefix , zpad , gflag ) ;
   } else if( tta2mni ){
     outset = THD_warp3D_tta2mni( inset , newggg , prefix,zpad,gflag ) ;
   } else if( mni2tta ){
     outset = THD_warp3D_mni2tta( inset , newggg , prefix,zpad,gflag ) ;
   }

   if( outset == NULL )
     ERROR_exit("THD_warp3D() fails for some reason!\n") ;

   /*-- polish up the new dataset info --*/
   /* if deobliquing, clear the oblique transformation matrix */
   /* really should update with new info, but clear for now */
   /* make invalid by setting lower right element to 0 */
#if 0
   if(oblique_flag==1) { 
     ZERO_MAT44(outset->daxes->ijk_to_dicom_real);
     outset->daxes->ijk_to_dicom_real.m[0][0] = 0.0;
   }
#endif

   /* lose obliquity if using 3dWarp for any transformation */
   /* recompute Tc (Cardinal transformation matrix for new grid output */
   THD_make_cardinal(outset);

   tross_Copy_History( inset , outset ) ;
   tross_Make_History( "3dWarp" , argc,argv , outset ) ;

   /*- 11 Mar 2004: mark output dataset as being in MNI space -*/

   if( tta2mni ){
/*     THD_set_string_atr( outset->dblk , "TLRC_SPACE" , "MNI-152" ) ;*/
     MCW_strncpy(outset->atlas_space, "MNI", THD_MAX_NAME);
   }


   if( mni2tta ){
     MCW_strncpy(outset->atlas_space, "TLRC", THD_MAX_NAME);
   }

   /*- 03 Aug 2005: save WARPDRIVE attributes, if present -*/

   if( atr_matfor != NULL ) {
      THD_insert_atr( outset->dblk, (ATR_any *)atr_matfor );
   } else if( use_matvec ){ /* put in the new matrix ZSS July 19 07*/
      UNLOAD_MAT(dicom_out2in.mm,matar[0],matar[1],matar[2],
                            matar[4],matar[5],matar[6],
                            matar[8],matar[9],matar[10] ) ;
       UNLOAD_FVEC3(dicom_out2in.vv,matar[3],matar[7],matar[11]) ;
       sprintf(anam,"WARPDRIVE_MATVEC_FOR_%06d",0) ;
       THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;

   }
   if( atr_matinv != NULL ) {
      THD_insert_atr( outset->dblk, (ATR_any *)atr_matinv );
   } else if( use_matvec ){ /* put in the new matrix ZSS July 19 07*/
       UNLOAD_MAT(dicom_in2out.mm,matar[0],matar[1],matar[2],
                            matar[4],matar[5],matar[6],
                            matar[8],matar[9],matar[10] ) ;
       UNLOAD_FVEC3(dicom_in2out.vv,matar[3],matar[7],matar[11]) ;
       sprintf(anam,"WARPDRIVE_MATVEC_INV_%06d",0) ;
       THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
   }
   /*- actually do the writositing -*/

   DSET_delete( inset ) ;
   DSET_write( outset ) ; if( verb ) WROTE_DSET(outset) ;
   exit(0) ;
}


#if 0
#define MAXNUM(a,b) ( (a) > (b) ? (a):(b))
#define MAX3(a,b,c) ( (MAXNUM(a,b)) > (MAXNUM(a,c)) ? (MAXNUM(a,b)):(MAXNUM(a,c)))
#define MINNUM(a,b) ( (a) < (b) ? (a):(b))
#define MIN3(a,b,c) ( (MINNUM(a,b)) < (MINNUM(a,c)) ? (MINNUM(a,b)):(MINNUM(a,c)))
/* compute angle of greatest obliquity given transformation matrix */
float compute_oblique_angle(mat44 ijk_to_dicom44)
{
   float dxtmp, dytmp, dztmp ;
   float xmax, ymax, zmax ;
   float fig_merit, ang_merit ;


   dxtmp = sqrt ( ijk_to_dicom44.m[0][0] * ijk_to_dicom44.m[0][0] +
                  ijk_to_dicom44.m[1][0] * ijk_to_dicom44.m[1][0] +
                  ijk_to_dicom44.m[2][0] * ijk_to_dicom44.m[2][0] ) ;

   xmax = MAX3(fabs(ijk_to_dicom44.m[0][0]),fabs(ijk_to_dicom44.m[1][0]),fabs(ijk_to_dicom44.m[2][0])) / dxtmp ;

   dytmp = sqrt ( ijk_to_dicom44.m[0][1] * ijk_to_dicom44.m[0][1] +
                  ijk_to_dicom44.m[1][1] * ijk_to_dicom44.m[1][1] +
                  ijk_to_dicom44.m[2][1] * ijk_to_dicom44.m[2][1] ) ;

   ymax = MAX3(fabs(ijk_to_dicom44.m[0][1]),
               fabs(ijk_to_dicom44.m[1][1]),
               fabs(ijk_to_dicom44.m[2][1])) / dytmp ;

   dztmp = sqrt ( ijk_to_dicom44.m[0][2] * ijk_to_dicom44.m[0][2] +
                  ijk_to_dicom44.m[1][2] * ijk_to_dicom44.m[1][2] +
                  ijk_to_dicom44.m[2][2] * ijk_to_dicom44.m[2][2] ) ;

   zmax = MAX3(fabs(ijk_to_dicom44.m[0][2]),
               fabs(ijk_to_dicom44.m[1][2]),
               fabs(ijk_to_dicom44.m[2][2])) / dztmp ;

   fig_merit = MIN3(xmax,ymax,zmax) ;
   ang_merit = acos (fig_merit) * 180.0 / 3.141592653 ;

   if (fabs(ang_merit) > .01) {
     INFO_message("%f degrees from plumb.\n",ang_merit ) ;
   }
   else 
      ang_merit = 0.0;
   return(ang_merit);
}
#endif


/* compute the transformation for deobliquing a dataset */
/* the IJK_TO_DICOM_REAL matrix must be stored in the dataset structure */
static void
Compute_Deoblique_Transformation(THD_3dim_dataset *dset, mat44 *Tw)
{
   THD_dmat33 tmat ;
   THD_dfvec3 tvec ;
   mat44 Tw_temp, Tc, inv_Tc, Tr;

   ENTRY("Compute_Deoblique_Transformation");
   if(!ISVALID_MAT44(dset->daxes->ijk_to_dicom_real)) {
      ERROR_exit("Oblique parent dataset doesn't have oblique "
                 "transformation attributes!?\n");
   }
   /* load oblique transformation matrix */
   Tr = dset->daxes->ijk_to_dicom_real;

   /* load DICOM (RAI) cardinal transformation matrix */
   THD_dicom_card_xform(dset, &tmat, &tvec); 
   LOAD_MAT44(Tc, tmat.mat[0][0], tmat.mat[0][1], tmat.mat[0][2], tvec.xyz[0],
                  tmat.mat[1][0], tmat.mat[1][1], tmat.mat[1][2], tvec.xyz[1],
		  tmat.mat[2][0], tmat.mat[2][1], tmat.mat[2][2], tvec.xyz[2]);
   inv_Tc = MAT44_INV(Tc);
   Tw_temp = MAT44_MUL(Tr, inv_Tc);
   *Tw = Tw_temp;

#ifdef DEBUG_ON
   DUMP_MAT44("Tr",Tr);
   DUMP_MAT44("Tc",Tc);
   DUMP_MAT44("Tw",Tw_temp);
#endif

   EXRETURN;
}
