00013 {
00014    int iarg=1 , nn , invert,nadd ;
00015    THD_dmat33 tmat , qmat , imat ;
00016    THD_dfvec3 tvec , qvec , ivec ;
00017    FILE * fp ;
00018    double dd[12] ;
00019    THD_dvecmat dvm ;
00020 
00021    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00022       printf("Usage: cat_matvec matvec_spec matvec_spec ...\n"
00023              "\n"
00024              "Catenates 3D rotation+shift matrix+vector transformations.\n"
00025              "Each matvec_spec is of the form\n"
00026              "\n"
00027              "  mfile [-opkey]\n"
00028              "\n"
00029              "'mfile' can take 2 forms:\n"
00030              "\n"
00031              "=== FORM 1 ===\n"
00032              "mfile is the name of an ASCII file with 12 numbers arranged\n"
00033              "in 3 lines:\n"
00034              "      u11 u12 u13 v1\n"
00035              "      u21 u22 u23 v2\n"
00036              "      u31 u32 u33 u3\n"
00037              "where each 'uij' and 'vi' is a number.  The 3x3 matrix [uij]\n"
00038              "is the matrix of the transform, and the 3-vector [vi] is the\n"
00039              "shift.  The transform is [xnew] = [uij]*[xold] + [vi].\n"
00040              "\n"
00041              "=== FORM 2 ===\n"
00042              "mfile is of the form 'dataset::attribute', where 'dataset'\n"
00043              "is the name of an AFNI dataset, and 'attribute' is the name\n"
00044              "of an attribute in the dataset's header that contains a\n"
00045              "matrix+vector (e.g., 'fred+orig::VOLREG_MATVEC_000000').\n"
00046              "\n"
00047              "=== COMPUTATIONS ===\n"
00048              "If [U] [v] are the matrix/vector for the first mfile, and\n"
00049              "   [A] [b] are the matrix/vector for the second mfile, then\n"
00050              "the catenated transformation is\n"
00051              "  matrix = [A][U]   vector = [A][v] + [b]\n"
00052              "That is, the second mfile transformation follows the first.\n"
00053              "\n"
00054              "The optional 'opkey' (operation key) following each mfile\n"
00055              "starts with a '-', and then is a set of letters telling how\n"
00056              "to treat the input.  The only opkey currently defined is\n"
00057              "\n"
00058              "  -I = invert the transformation:\n"
00059              "                     -1              -1\n"
00060              "       [xold] = [uij]  [xnew] - [uij]  [vi]\n"
00061              "\n"
00062              "The transformation resulting by catenating the transformations\n"
00063              "is written to stdout in the sam ASCII file format.  This can be\n"
00064              "used as input to 3drotate -matvec_dicom (provided [uij] is a\n"
00065              "proper orthogonal matrix).\n"
00066              "\n"
00067              "N.B.: If exactly 9 numbers can be read from an mfile, then those\n"
00068              "      values form the matrix, and the vector is set to zero.\n"
00069            ) ;
00070       exit(0) ;
00071    }
00072 
00073    /* initialize identity transformation into tmat,tvec */
00074 
00075    LOAD_DIAG_DMAT(tmat,1.0,1.0,1.0) ;
00076    LOAD_DFVEC3(tvec,0.0,0.0,0.0) ;
00077 
00078    /* loop and read arguments, process them */
00079 
00080    while( iarg < argc ){
00081 
00082       nadd = 1 ; invert = 0 ;
00083       if( iarg+1 < argc && strcmp(argv[iarg+1],"-I") == 0 ){
00084          invert = 1 ; nadd = 2 ;
00085       }
00086       dvm = THD_read_dvecmat( argv[iarg] , invert ) ;
00087       qmat = dvm.mm ; qvec = dvm.vv ;
00088 
00089       if( SIZE_DMAT(qmat) == 0.0 ){
00090          fprintf(stderr,"*** ERROR: can't read mfile %s\n",argv[iarg]) ;
00091          exit(1) ;
00092       }
00093 
00094       iarg += nadd ;
00095 
00096       /* multiply into accumulating transformation */
00097 
00098       imat = DMAT_MUL( qmat , tmat ) ;
00099       ivec = DMATVEC( qmat , tvec ) ;
00100       tvec = ADD_DFVEC3( qvec , ivec ) ;
00101       tmat = imat ;
00102    }
00103 
00104    /* write results to stdout */
00105 
00106    printf("%13.6g %13.6g %13.6g %13.6g\n"
00107           "%13.6g %13.6g %13.6g %13.6g\n"
00108           "%13.6g %13.6g %13.6g %13.6g\n" ,
00109      tmat.mat[0][0] , tmat.mat[0][1] , tmat.mat[0][2] , tvec.xyz[0] ,
00110      tmat.mat[1][0] , tmat.mat[1][1] , tmat.mat[1][2] , tvec.xyz[1] ,
00111      tmat.mat[2][0] , tmat.mat[2][1] , tmat.mat[2][2] , tvec.xyz[2]  ) ;
00112 
00113    exit(0) ;
00114 }