00001 #include "mrilib.h"
00002 
00003 /*-------------------------------------------------------------------
00004   08 Feb 2001: read a matrix+vector from a file - RWCox
00005   14 Feb 2001: modified to add "-rotate ..." mode
00006 ---------------------------------------------------------------------*/
00007 
00008 THD_dvecmat THD_read_dvecmat( char * fname , int invert )
00009 {
00010    THD_dvecmat dvm ;
00011    THD_dmat33 tmat ;
00012    THD_dfvec3 tvec ;
00013    FILE * fp ;
00014    double dd[12] ;
00015    int  nn ;
00016 
00017 ENTRY("THD_read_dvecmat") ;
00018 
00019    LOAD_DIAG_DMAT(tmat,0.0,0.0,0.0) ;  /* initialize to all 0s */
00020    LOAD_DFVEC3(tvec,0.0,0.0,0.0) ;
00021 
00022    if( fname == NULL || fname[0] == '\0' ){
00023       dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
00024    }
00025 
00026    /*-- filename is "dataset::attribute" --*/
00027 
00028    if( strstr(fname,"::") != NULL ){  /*== read from dataset header ==*/
00029       char * dname = strdup(fname) ;
00030       char * cc = strstr(dname,"::") ;
00031       THD_3dim_dataset * dset ;
00032       ATR_float * atr ;
00033 
00034       *cc = '\0' ; cc += 2 ;  /* dname = dataset name ; cc = attribute name */
00035 
00036       dset = THD_open_one_dataset( dname ) ;
00037       if( !ISVALID_DSET(dset) ){
00038          fprintf(stderr,
00039                  "*** THD_read_dvecmat: can't open dataset %s\n",
00040                  dname) ;
00041          dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
00042       }
00043 
00044       atr = THD_find_float_atr( dset->dblk , cc ) ;
00045       if( atr == NULL ){
00046          fprintf(stderr,
00047                  "*** THD_read_dvecmat: can't find attribute %s in dataset %s\n",
00048                  cc,dname) ;
00049          dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
00050       }
00051 
00052       switch( atr->nfl ){
00053 
00054          default:
00055             fprintf(stderr,
00056                     "*** THD_read_dvecmat: can't read matrix+vector from dataset %s, attribute %s\n",
00057                     dname,cc) ;
00058             dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
00059          break ; /* unreachable */
00060 
00061          case 12: LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
00062                                  atr->fl[4],atr->fl[5],atr->fl[6],
00063                                  atr->fl[8],atr->fl[9],atr->fl[10] ) ;
00064                   LOAD_DFVEC3(tvec,atr->fl[3],atr->fl[7],atr->fl[11]) ;
00065          break ;
00066 
00067          case 9:  LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
00068                                  atr->fl[3],atr->fl[4],atr->fl[5],
00069                                  atr->fl[6],atr->fl[7],atr->fl[8] ) ;
00070                   LOAD_DFVEC3(tvec,0,0,0) ;
00071          break ;
00072       }
00073 
00074       free(dname) ; DSET_delete(dset) ;
00075 
00076    /*-- 14 Feb 2001: filename is "-rotate a b c -[ab]shift x y z" string --*/
00077 
00078    } else if( strstr(fname,"-rotate") != NULL ){  /*== compute directly ==*/
00079 
00080       dvm = THD_rotcom_to_matvec( NULL , fname ) ;  /* thd_rotangles.c */
00081       tvec = dvm.vv ; tmat = dvm.mm ;
00082 
00083    /*-- just a normal filename --*/
00084 
00085    } else {  /*== read numbers from file ==*/
00086 
00087       fp = fopen( fname , "r" ) ;
00088       if( fp == NULL ){
00089          fprintf(stderr,
00090                  "*** THD_read_dvecmat: can't open file %s\n",
00091                  fname) ;
00092          dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
00093       }
00094 
00095       nn = fscanf(fp,"%lf %lf %lf %lf"      /* try to get 12 numbers */
00096                      "%lf %lf %lf %lf"
00097                      "%lf %lf %lf %lf" ,
00098                   dd+0,dd+1,dd+2 ,dd+3,
00099                   dd+4,dd+5,dd+6 ,dd+7,
00100                   dd+8,dd+9,dd+10,dd+11 ) ;
00101 
00102       switch( nn ){  /* how many did we actually read? */
00103 
00104          default:
00105             fprintf(stderr,
00106                     "*** THD_read_dvecmat: can't read matrix+vector from file %s\n",
00107                     fname) ;
00108             dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
00109          break ; /* unreachable */
00110 
00111          case 12: LOAD_DMAT(tmat,dd[0],dd[1],dd[2],    /* 12 ==> matrix+vector */
00112                                  dd[4],dd[5],dd[6],
00113                                  dd[8],dd[9],dd[10] ) ;
00114                   LOAD_DFVEC3(tvec,dd[3],dd[7],dd[11]) ;
00115          break ;
00116 
00117          case 9:  LOAD_DMAT(tmat,dd[0],dd[1],dd[2],    /*  9 ==> matrix only */
00118                                  dd[3],dd[4],dd[5],
00119                                  dd[6],dd[7],dd[8] ) ;
00120                   LOAD_DFVEC3(tvec,0,0,0) ;
00121          break ;
00122       }
00123    }
00124 
00125    /*-- invert the transformation we just read?            --*/
00126    /*-- [y] = [R][x] + [v]       is transformation, so     --*/
00127    /*-- [x] = inv[R] - inv[R][v] is inverse transformation --*/
00128 
00129    if( invert ){
00130       THD_dmat33 imat ; THD_dfvec3 ivec ;
00131       imat = DMAT_INV(tmat) ;             /* matrix inverse */
00132       ivec = DMATVEC(imat,tvec) ;         /* multiply inverse into vector */
00133       tmat = imat ;
00134       tvec = ivec ; NEGATE_DFVEC3(tvec) ; /* negate vector */
00135    }
00136 
00137    /*-- store results and get outta here, dude! --*/
00138 
00139    dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
00140 }