#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 ***/
static THD_fvec3 AFNI_forward_warp_vector ( THD_warp * , THD_fvec3 ) ;
static THD_fvec3 AFNI_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"
) ;
exit(0) ;
}
/*--------------------------------------------------------------------------*/
static void errex( char * str )
{
fprintf(stderr,"** FATAL ERROR: %s\n",str) ; exit(1) ;
}
/*--------------------------------------------------------------------------*/
#define NBUF 1024
#define isnumeric(c) (isdigit(c) || c == '-' || c == '+' || c == '.')
#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 ;
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 );
oset = THD_open_dataset(lbuf) ;
if( !ISVALID_DSET(oset) ){
char str[NBUF] ;
sprintf(str,"-apar: Can't open dataset %s",lbuf) ;
errex(str) ;
}
warp = aset->warp ;
iarg++ ; continue ;
}
/* -matvec */
if( strcmp(argv[iarg],"-matvec") == 0 ){
THD_dvecmat dvm ;
THD_linear_mapping 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 = 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 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 = 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 = 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 = 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 = AFNI_backward_warp_vector( warp , vin ) ;
else vout = AFNI_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 = 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
--------------------------------------------------------------------------*/
static THD_fvec3 AFNI_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
--------------------------------------------------------------------------*/
static THD_fvec3 AFNI_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
syntax highlighted by Code2HTML, v. 0.9.1