#include "mrilib.h"
/*-------------------------------------------------------------------
08 Feb 2001: read a matrix+vector from a file - RWCox
14 Feb 2001: modified to add "-rotate ..." mode
10 Aug 2005: add MATRIX(...)
---------------------------------------------------------------------*/
THD_dvecmat THD_read_dvecmat( char *fname , int invert )
{
THD_dvecmat dvm ;
THD_dmat33 tmat ;
THD_dfvec3 tvec ;
FILE *fp ;
double dd[12] ;
int nn , loaded=0 ;
ENTRY("THD_read_dvecmat") ;
LOAD_DIAG_DMAT(tmat,0.0,0.0,0.0) ; /* initialize to all 0s */
LOAD_DFVEC3(tvec,0.0,0.0,0.0) ;
if( fname == NULL || fname[0] == '\0' ){
dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
}
/*-- filename is "dataset::attribute" --*/
if( strstr(fname,"::") != NULL ){ /*== read from dataset header ==*/
char *dname = strdup(fname) ;
char *cc = strstr(dname,"::") ;
char *ss ; int iss=0 ;
THD_3dim_dataset *dset ;
ATR_float *atr, *atrc ;
int incl;
THD_dfvec3 tvec_co, tvec_cb;
*cc = '\0' ; cc += 2 ; /* dname = dataset name ; cc = attribute name */
dset = THD_open_one_dataset( dname ) ;
if( !ISVALID_DSET(dset) ){
ERROR_message("THD_read_dvecmat: can't open dataset %s\n",dname) ;
dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
}
if (strstr(cc,"IJK_TO_CARD_DICOM")) { /* return the matrix that turns i,j,k to mm RAI ZSS May 07*/
THD_dicom_card_xform(dset, &tmat, &tvec);
} else {
ss = strstr(cc,"[") ;
if( ss != NULL ){
*ss = '\0'; ss++; iss = (int)strtol(ss,NULL,10); if( iss < 0 ) iss = 0;
}
atr = THD_find_float_atr( dset->dblk , cc ) ;
if( atr == NULL ){
ERROR_message("THD_read_dvecmat: can't find attribute %s in dataset %s\n",
cc,dname) ;
dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
}
if( strcmp(cc,"WARP_DATA") == 0 ){ /* 10 Aug 2005 */
LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
atr->fl[3],atr->fl[4],atr->fl[5],
atr->fl[6],atr->fl[7],atr->fl[8] ) ;
LOAD_DFVEC3(tvec,-atr->fl[18],-atr->fl[19],-atr->fl[20]) ;
loaded = 1 ;
}
if( !loaded ){
switch( atr->nfl ){
default:
ERROR_message(
"THD_read_dvecmat: can't read matrix+vector from %s::%s\n",
dname,cc) ;
dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
break ; /* unreachable */
case 12: LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
atr->fl[4],atr->fl[5],atr->fl[6],
atr->fl[8],atr->fl[9],atr->fl[10] ) ;
LOAD_DFVEC3(tvec,atr->fl[3],atr->fl[7],atr->fl[11]) ;
break ;
case 9: LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
atr->fl[3],atr->fl[4],atr->fl[5],
atr->fl[6],atr->fl[7],atr->fl[8] ) ;
LOAD_DFVEC3(tvec,0,0,0) ;
break ;
}
}
{
/* ZSS: Dec 27, the year of our Lord when the War on Christmas was raging */
/* Need to include VOLREG_CENTER_OLD, VOLREG_CENTER_BASE,
ROTATE_CENTER_OLD, ROTATE_CENTER_BASE, if applicable. */
/* Do we have a ROTATE or VOLREG attribute ?*/
incl = 0;
if (strstr(cc,"VOLREG_MATVEC")) {
atrc = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_OLD");
if ( atrc ) {
LOAD_DFVEC3(tvec_co,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
incl = 1;
} else {
LOAD_DFVEC3(tvec_co,0,0,0) ;
}
atrc = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_BASE");
if ( atrc ) {
LOAD_DFVEC3(tvec_cb,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
incl = 1;
} else {
LOAD_DFVEC3(tvec_cb,0,0,0) ;
}
if (incl == 1) INFO_message("THD_read_dvecmat:\n"
" Including VOLREG_CENTER_BASE and VOLREG_CENTER_OLD\n"
" attributes in final transform\n");
else INFO_message("THD_read_dvecmat:\n"
" No VOLREG_CENTER_BASE or VOLREG_CENTER_OLD\n"
" attributes found with VOLREG_MATVEC\n");
} else if (strstr(cc,"ROTATE_MATVEC")) {
atrc = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_OLD");
if ( atrc ) {
LOAD_DFVEC3(tvec_co,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
incl = 1;
} else {
LOAD_DFVEC3(tvec_co,0,0,0) ;
}
atrc = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_BASE");
if ( atrc ) {
LOAD_DFVEC3(tvec_cb,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
incl = 1;
} else {
LOAD_DFVEC3(tvec_cb,0,0,0) ;
}
if (incl == 1) INFO_message("THD_read_dvecmat:\n"
" Including ROTATE_CENTER_BASE and ROTATE_CENTER_OLD\n"
" attributes in final transform\n");
else INFO_message("THD_read_dvecmat:\n"
" No ROTATE_CENTER_BASE or ROTATE_CENTER_OLD\n"
" attributes found with ROTATE_MATVEC\n");
}
if (incl == 1) {
tvec_co = DMATVEC(tmat, tvec_co);
NEGATE_DFVEC3(tvec_co);
tvec.xyz[0] += tvec_cb.xyz[0] + tvec_co.xyz[0];
tvec.xyz[1] += tvec_cb.xyz[1] + tvec_co.xyz[1];
tvec.xyz[2] += tvec_cb.xyz[2] + tvec_co.xyz[2];
}
}
}
free(dname) ; DSET_delete(dset) ;
/*-- 14 Feb 2001: filename is "-rotate a b c -[ab]shift x y z" string --*/
} else if( strstr(fname,"-rotate") != NULL ){ /*== compute directly ==*/
dvm = THD_rotcom_to_matvec( NULL , fname ) ; /* thd_rotangles.c */
tvec = dvm.vv ; tmat = dvm.mm ;
/*-- MATRIX(...) --*/
} else if( strncmp(fname,"MATRIX(",7) == 0 ){
float matar[12] ; int nn ;
nn = sscanf(fname,"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_message("badly formatted: %s",fname) ;
} else {
LOAD_DMAT(tmat,matar[0],matar[1],matar[2],
matar[4],matar[5],matar[6],
matar[8],matar[9],matar[10] ) ;
LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
}
/*-- just a normal filename --*/
} else { /*== read numbers from file ==*/
fp = fopen( fname , "r" ) ;
if( fp == NULL ){
ERROR_message(
"THD_read_dvecmat: can't open file %s\n",
fname) ;
dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
}
nn = fscanf(fp,"%lf %lf %lf %lf" /* try to get 12 numbers */
"%lf %lf %lf %lf"
"%lf %lf %lf %lf" ,
dd+0,dd+1,dd+2 ,dd+3,
dd+4,dd+5,dd+6 ,dd+7,
dd+8,dd+9,dd+10,dd+11 ) ;
switch( nn ){ /* how many did we actually read? */
default:
ERROR_message(
"THD_read_dvecmat: can't read matrix+vector from file %s\n",
fname) ;
dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
break ; /* unreachable */
case 12: LOAD_DMAT(tmat,dd[0],dd[1],dd[2], /* 12 ==> matrix+vector */
dd[4],dd[5],dd[6],
dd[8],dd[9],dd[10] ) ;
LOAD_DFVEC3(tvec,dd[3],dd[7],dd[11]) ;
break ;
case 9: LOAD_DMAT(tmat,dd[0],dd[1],dd[2], /* 9 ==> matrix only */
dd[3],dd[4],dd[5],
dd[6],dd[7],dd[8] ) ;
LOAD_DFVEC3(tvec,0,0,0) ;
break ;
}
}
/*-- invert the transformation we just read? --*/
/*-- [y] = [R][x] + [v] is transformation, so --*/
/*-- [x] = inv[R] - inv[R][v] is inverse transformation --*/
if( invert ){
THD_dmat33 imat ; THD_dfvec3 ivec ;
imat = DMAT_INV(tmat) ; /* matrix inverse */
ivec = DMATVEC(imat,tvec) ; /* multiply inverse into vector */
tmat = imat ;
tvec = ivec ; NEGATE_DFVEC3(tvec) ; /* negate vector */
}
/*-- store results and get outta here, dude! --*/
dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
}
syntax highlighted by Code2HTML, v. 0.9.1