00505 {
00506    THD_fvec3 new_fv ;
00507 
00508    if( warp == NULL ) return old_fv ;
00509 
00510    switch( warp->type ){
00511 
00512       default: new_fv = old_fv ; break ;
00513 
00514       case WARP_TALAIRACH_12_TYPE:{
00515          THD_linear_mapping map ;
00516          int iw ;
00517 
00518          /* test if input is in bot..top of each defined map */
00519 
00520          for( iw=0 ; iw < 12 ; iw++ ){
00521             map = warp->tal_12.warp[iw] ;
00522 
00523             if( old_fv.xyz[0] >= map.bot.xyz[0] &&
00524                 old_fv.xyz[1] >= map.bot.xyz[1] &&
00525                 old_fv.xyz[2] >= map.bot.xyz[2] &&
00526                 old_fv.xyz[0] <= map.top.xyz[0] &&
00527                 old_fv.xyz[1] <= map.top.xyz[1] &&
00528                 old_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
00529          }
00530          new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00531       }
00532       break ;
00533 
00534       case WARP_AFFINE_TYPE:{
00535          THD_linear_mapping map = warp->rig_bod.warp ;
00536          new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00537       }
00538       break ;
00539 
00540    }
00541    return new_fv ;
00542 }

00460 {
00461    THD_fvec3 new_fv ;
00462 
00463    if( warp == NULL ) return old_fv ;
00464 
00465    switch( warp->type ){
00466 
00467       default: new_fv = old_fv ; break ;
00468 
00469       case WARP_TALAIRACH_12_TYPE:{
00470          THD_linear_mapping map ;
00471          int iw ;
00472 
00473          /* forward transform each possible case,
00474             and test if result is in bot..top of defined map */
00475 
00476          for( iw=0 ; iw < 12 ; iw++ ){
00477             map    = warp->tal_12.warp[iw] ;
00478             new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00479 
00480             if( new_fv.xyz[0] >= map.bot.xyz[0] &&
00481                 new_fv.xyz[1] >= map.bot.xyz[1] &&
00482                 new_fv.xyz[2] >= map.bot.xyz[2] &&
00483                 new_fv.xyz[0] <= map.top.xyz[0] &&
00484                 new_fv.xyz[1] <= map.top.xyz[1] &&
00485                 new_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
00486          }
00487       }
00488       break ;
00489 
00490       case WARP_AFFINE_TYPE:{
00491          THD_linear_mapping map = warp->rig_bod.warp ;
00492          new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00493       }
00494       break ;
00495 
00496    }
00497    return new_fv ;
00498 }

00151 {
00152    fprintf(stderr,"** FATAL ERROR: %s\n",str) ; exit(1) ;
00153 }

00164 {
00165    int iarg=1 ;
00166    FILE *fin=stdin , *fout=stdout ;
00167    int backward=0 , force=0 ;
00168    THD_warp *warp=NULL ;
00169    char lbuf[NBUF] , *cpt ;
00170    int itype=AFNI_1D , numv=0 , numc=0 ;
00171    THD_fvec3 vin , vout ;
00172    float xx,yy,zz ;
00173    int nn , ii , good=0 ;
00174    THD_3dim_dataset *aset=NULL , *oset=NULL ;
00175 
00176    if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Syntax() ;
00177 
00178    /*-- process command line arguments --*/
00179 
00180    while( iarg < argc ){
00181 
00182       /* -input */
00183 
00184       if( strcmp(argv[iarg],"-input") == 0 ){
00185          if( ++iarg >= argc )
00186             errex("-input: Need argument after -input") ;
00187          if( fin != stdin )
00188             errex("-input: Can't have two -input options") ;
00189          if( strcmp(argv[iarg],"-") != 0 ){
00190             fin = fopen( argv[iarg] , "r" ) ;
00191             if( fin == NULL )
00192                errex("-input: Can't open input file") ;
00193          }
00194          iarg++ ; continue ;
00195       }
00196 
00197       /* -output */
00198 
00199       if( strcmp(argv[iarg],"-output") == 0 ){
00200          if( ++iarg >= argc )
00201             errex("-output: Need argument after -output") ;
00202          if( fout != stdout )
00203             errex("-output: Can't have two -output options") ;
00204          if( strcmp(argv[iarg],"-") != 0 ){
00205             if( !THD_filename_ok(argv[iarg]) )
00206                errex("-output: Output filename is illegal") ;
00207             if( !force && THD_is_file(argv[iarg]) )
00208                errex("-output: Output file already exists") ;
00209             fout = fopen( argv[iarg] , "w" ) ;
00210             if( fout == NULL )
00211                errex("-output: Can't open output file") ;
00212          }
00213          iarg++ ; continue ;
00214       }
00215 
00216       /* -force */
00217 
00218       if( strcmp(argv[iarg],"-force") == 0 ){
00219          force = 1 ;
00220          iarg++ ; continue ;
00221       }
00222 
00223       /* -forward */
00224 
00225       if( strcmp(argv[iarg],"-forward") == 0 ){
00226          backward = 0 ;
00227          iarg++ ; continue ;
00228       }
00229 
00230       /* -backward */
00231 
00232       if( strcmp(argv[iarg],"-backward") == 0 ){
00233          backward = 1 ;
00234          iarg++ ; continue ;
00235       }
00236 
00237       /* -apar */
00238 
00239       if( strcmp(argv[iarg],"-apar") == 0 ){
00240          if( ++iarg >= argc )
00241             errex("-apar: Need argument after -apar") ;
00242          if( warp != NULL )
00243             errex("-apar: Can't specify transformation twice") ;
00244 
00245          /* open dataset with warp */
00246 
00247          aset = THD_open_dataset( argv[iarg] ) ;
00248          if( !ISVALID_DSET(aset) ){
00249             sprintf(lbuf,"-apar: Can't open dataset %s\n",argv[iarg]) ;
00250             errex(lbuf) ;
00251          }
00252          if( aset->warp == NULL ){
00253             sprintf(lbuf,"-apar: Dataset %s does not contain warp",argv[iarg]) ;
00254             errex(lbuf) ;
00255          }
00256          if( aset->view_type == VIEW_ORIGINAL_TYPE ){
00257             sprintf(lbuf,"-apar: Dataset %s is in the +orig view",argv[iarg]) ;
00258             errex(lbuf) ;
00259          }
00260 
00261          /* open +orig version of this dataset */
00262 
00263          sprintf(lbuf,"%s%s+orig.HEAD", aset->dblk->diskptr->directory_name,
00264                                         aset->dblk->diskptr->prefix         );
00265 
00266          oset = THD_open_dataset(lbuf) ;
00267          if( !ISVALID_DSET(oset) ){
00268             char str[NBUF] ;
00269             sprintf(str,"-apar: Can't open dataset %s",lbuf) ;
00270             errex(str) ;
00271          }
00272 
00273          warp = aset->warp ;
00274          iarg++ ; continue ;
00275       }
00276 
00277       /* -matvec */
00278 
00279       if( strcmp(argv[iarg],"-matvec") == 0 ){
00280          THD_dvecmat dvm ;
00281          THD_linear_mapping lmap ;
00282 
00283          if( ++iarg >= argc )
00284             errex("-matvec: Need argument after -matvec") ;
00285          if( warp != NULL )
00286             errex("-matvec: Can't specify transformation twice") ;
00287          dvm = THD_read_dvecmat( argv[iarg] , 0 ) ;
00288          if( SIZE_DMAT(dvm.mm) == 0.0 )
00289             errex("-matvec: Can't read transformation") ;
00290 
00291          /* manufacture an AFNI linear map */
00292 
00293          lmap.type = MAPPING_LINEAR_TYPE ;
00294          lmap.mfor.mat[0][0] = dvm.mm.mat[0][0] ;  /* copy matrix in */
00295          lmap.mfor.mat[0][1] = dvm.mm.mat[0][1] ;
00296          lmap.mfor.mat[0][2] = dvm.mm.mat[0][2] ;
00297          lmap.mfor.mat[1][0] = dvm.mm.mat[1][0] ;
00298          lmap.mfor.mat[1][1] = dvm.mm.mat[1][1] ;
00299          lmap.mfor.mat[1][2] = dvm.mm.mat[1][2] ;
00300          lmap.mfor.mat[2][0] = dvm.mm.mat[2][0] ;
00301          lmap.mfor.mat[2][1] = dvm.mm.mat[2][1] ;
00302          lmap.mfor.mat[2][2] = dvm.mm.mat[2][2] ;
00303          lmap.bvec.xyz[0]    = -dvm.vv.xyz[0] ;    /* copy vector in */
00304          lmap.bvec.xyz[1]    = -dvm.vv.xyz[1] ;
00305          lmap.bvec.xyz[2]    = -dvm.vv.xyz[2] ;
00306          LOAD_INVERSE_LMAP(lmap) ;                 /* make inverse map */
00307 
00308          /* manufacture an AFNI warp */
00309 
00310          warp = (THD_warp *) calloc( sizeof(THD_warp) , 1 ) ;
00311          warp->type = WARP_AFFINE_TYPE ;
00312          warp->rig_bod.warp = lmap ;
00313 
00314          fprintf(stderr,"++ Using matrix-vector transformation below:\n"
00315                         "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n"
00316                         "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n"
00317                         "   [ %9.5f %9.5f %9.5f ]   [ %9.5f ]\n" ,
00318            dvm.mm.mat[0][0] , dvm.mm.mat[0][1] , dvm.mm.mat[0][2] ,
00319            dvm.mm.mat[1][0] , dvm.mm.mat[1][1] , dvm.mm.mat[1][2] ,
00320            dvm.mm.mat[2][0] , dvm.mm.mat[2][1] , dvm.mm.mat[2][2] ,
00321            dvm.vv.xyz[0]    , dvm.vv.xyz[1]    , dvm.vv.xyz[2]     ) ;
00322 
00323          iarg++ ; continue ;
00324       }
00325 
00326       /* ERROR */
00327 
00328       { char *str = AFMALL(char, strlen(argv[iarg])+32) ;
00329         sprintf(str,"Unknown option: %s",argv[iarg]) ;
00330         errex(str) ;
00331       }
00332 
00333    } /* end of loop over command line args */
00334 
00335    /*-- check if we are ready --*/
00336 
00337    if( warp == NULL )
00338       errex("Transformation not specified") ;
00339 
00340    /*-- Read 1st line of input file to determine what to do with it --*/
00341 
00342    cpt = fgets( lbuf , NBUF , fin ) ;
00343    if( cpt == NULL )
00344       errex("Couldn't read 1st line from input file") ;
00345 
00346    /*-- check if this is a SureFit file --*/
00347 
00348    if( strstr(cpt,"BeginHeader") != NULL ) itype = SUREFIT ;
00349    else                                    itype = AFNI_1D ;
00350 
00351    /*-- if SureFit, echo all header stuff to the output --*/
00352 
00353    if( itype == SUREFIT ){
00354       int numh=1 ;
00355       fprintf(fout,"%s",lbuf) ;    /* echo 1st line */
00356       while(1){                    /* read lines until EndHeader is found */
00357 
00358          cpt = fgets( lbuf , NBUF , fin ) ;
00359          if( cpt == NULL )
00360             errex("Input ended before EndHeader was found") ;
00361 
00362          fprintf(fout,"%s",lbuf) ; /* echo line to output */
00363          numh++ ;
00364 
00365          if( strstr(lbuf,"EndHeader") != NULL ){  /* do next line, too */
00366             cpt = fgets( lbuf , NBUF , fin ) ;
00367             if( cpt == NULL )
00368                errex("Input ended just after EndHeader") ;
00369             fprintf(fout,"%s",lbuf) ;
00370             fprintf(stderr,"++ Wrote %d SureFit header lines to output\n",numh);
00371             break ;                /* end of loop */
00372          }
00373       }
00374 
00375       cpt = fgets( lbuf , NBUF , fin ) ;  /* get next line, with 1st vector */
00376       if( cpt == NULL )
00377          errex("Input ended just after Node count") ;
00378 
00379    } /* end of SureFit header echo */
00380 
00381    /*-- From this point on, process each line with 1 vector. --*/
00382    /*-- At the start of the loop, 1 line's data is in lbuf.  --*/
00383 
00384    do{
00385 
00386       switch( itype ){
00387          case SUREFIT:
00388             ii = sscanf(lbuf,"%d%f%f%f",&nn,&xx,&yy,&zz) ;
00389             good = (ii == 4) ;
00390          break ;
00391 
00392          case AFNI_1D:
00393             ii = sscanf(lbuf,"%f%f%f",&xx,&yy,&zz) ;
00394             good = (ii == 3) ;
00395          break ;
00396       }
00397 
00398       if( !good ){  /* just echo line to output */
00399 
00400          fprintf(fout,"%s",lbuf) ;
00401          numc++ ;
00402 
00403       } else {      /* transform a vector!!! */
00404 
00405          LOAD_FVEC3(vin,xx,yy,zz) ;
00406 #ifdef DEBUG
00407 fprintf(stderr,"\n") ;
00408 DUMP_FVEC3("vin              ",vin) ;
00409 #endif
00410 
00411          if( itype == SUREFIT && aset != NULL ){
00412             if( backward ) vin = THD_surefit_to_dicomm( aset , vin ) ;
00413             else           vin = THD_surefit_to_dicomm( oset , vin ) ;
00414 #ifdef DEBUG
00415 DUMP_FVEC3("surefit_to_dicomm",vin) ;
00416 #endif
00417          }
00418 
00419          if( backward ) vout = AFNI_backward_warp_vector( warp , vin ) ;
00420          else           vout = AFNI_forward_warp_vector ( warp , vin ) ;
00421 
00422 #ifdef DEBUG
00423 DUMP_FVEC3("vout             ",vout) ;
00424 #endif
00425 
00426          if( itype == SUREFIT && aset != NULL ){
00427             if( backward ) vout = THD_dicomm_to_surefit( oset , vout ) ;
00428             else           vout = THD_dicomm_to_surefit( aset , vout ) ;
00429 #ifdef DEBUG
00430 DUMP_FVEC3("dicomm_to_surefit",vout) ;
00431 #endif
00432          }
00433 
00434          if( itype == SUREFIT ) fprintf(fout,"%d ",nn) ;
00435          fprintf(fout,"%f %f %f\n",vout.xyz[0],vout.xyz[1],vout.xyz[2]) ;
00436          numv++ ;
00437       }
00438 
00439       cpt = fgets( lbuf , NBUF , fin ) ;  /* get next line */
00440 
00441    } while( cpt != NULL ) ;  /* loop until no data can be read */
00442 
00443    if( numc > 0 )
00444       fprintf(stderr,"++ Wrote %d vectors;  %d comments\n",numv,numc) ;
00445    else
00446       fprintf(stderr,"++ Wrote %d vectors\n",numv) ;
00447 
00448    exit(0) ;
00449 }

00026 {
00027    printf(
00028     "Usage: Vecwarp [options]\n"
00029     "Transforms (warps) a list of 3-vectors into another list of 3-vectors\n"
00030     "according to the options.  Error messages, warnings, and informational\n"
00031     "messages are written to stderr.  If a fatal error occurs, the program\n"
00032     "exits with status 1; otherwise, it exits with status 0.\n"
00033     "\n"
00034     "OPTIONS:\n"
00035     " -apar aaa   = Use the AFNI dataset 'aaa' as the source of the\n"
00036     "               transformation; this dataset must be in +acpc\n"
00037     "               or +tlrc coordinates, and must contain the\n"
00038     "               attributes WARP_TYPE and WARP_DATA which describe\n"
00039     "               the forward transformation from +orig coordinates\n"
00040     "               to the 'aaa' coordinate system.\n"
00041     "             N.B.: The +orig version of this dataset must also be\n"
00042     "                   readable, since it is also needed when translating\n"
00043     "                   vectors between SureFit and AFNI coordinates.\n"
00044     "                   Only the .HEAD files are actually used.\n"
00045     "\n"
00046     " -matvec mmm = Read an affine transformation matrix-vector from file\n"
00047     "               'mmm', which must be in the format\n"
00048     "                   u11 u12 u13 v1\n"
00049     "                   u21 u22 u23 v2\n"
00050     "                   u31 u32 u33 v3\n"
00051     "               where each 'uij' and 'vi' is a number.  The forward\n"
00052     "               transformation is defined as\n"
00053     "                   [ xout ]   [ u11 u12 u13 ] [ xin ]   [ v1 ]\n"
00054     "                   [ yout ] = [ u21 u22 u23 ] [ yin ] + [ v2 ]\n"
00055     "                   [ zout ]   [ u31 u32 u33 ] [ zin ]   [ v3 ]\n"
00056     "\n"
00057     " Exactly one of -apar or -matvec must be used to specify the\n"
00058     " transformation.\n"
00059     "\n"
00060     " -forward    = -forward means to apply the forward transformation;\n"
00061     "   *OR*        -backward means to apply the backward transformation\n"
00062     " -backward     * For example, if the transformation is specified by\n"
00063     "                  '-apar fred+tlrc', then the forward transformation\n"
00064     "                  is from +orig to +tlrc coordinates, and the backward\n"
00065     "                  transformation is from +tlrc to +orig coordinates.\n"
00066     "               * If the transformation is specified by -matvec, then\n"
00067     "                  the matrix-vector read in defines the forward\n"
00068     "                  transform as above, and the backward transformation\n"
00069     "                  is defined as the inverse.\n"
00070     "               * If neither -forward nor -backward is given, then\n"
00071     "                  -forward is the default.\n"
00072     "\n"
00073     " -input iii  = Read input 3-vectors from file 'iii' (from stdin if\n"
00074     "               'iii' is '-' or the -input option is missing).  Input\n"
00075     "               data may be in one of the following ASCII formats:\n"
00076     "\n"
00077     "               * SureFit .coord files:\n"
00078     "                   BeginHeader\n"
00079     "                   lines of text ...\n"
00080     "                   EndHeader\n"
00081     "                   count\n"
00082     "                   int x y z\n"
00083     "                   int x y z\n"
00084     "                   et cetera...\n"
00085     "                 In this case, everything up to and including the\n"
00086     "                 count is simply passed through to the output.  Each\n"
00087     "                 (x,y,z) triple is transformed, and output with the\n"
00088     "                 int label that precedes it.  Lines that cannot be\n"
00089     "                 scanned as 1 int and 3 floats are treated as comments\n"
00090     "                 and are passed to through to the output unchanged.\n"
00091     "               N.B.: SureFit coordinates are\n"
00092     "                   x = distance Right    of Left-most      dataset corner\n"
00093     "                   y = distance Anterior to Posterior-most dataset corner\n"
00094     "                   z = distance Superior to Inferior-most  dataset corner\n"
00095     "                 For example, if the transformation is specified by\n"
00096     "                   -forward -apar fred+tlrc\n"
00097     "                 then the input (x,y,z) are relative to fred+orig and the\n"
00098     "                 output (x,y,z) are relative to fred+tlrc.  If instead\n"
00099     "                   -backward -apar fred+tlrc\n"
00100     "                 is used, then the input (x,y,z) are relative to fred+tlrc\n"
00101     "                 and the output (x,y,z) are relative to fred+orig.\n"
00102     "                 For this to work properly, not only fred+tlrc must be\n"
00103     "                 readable by Vecwarp, but fred+orig must be as well.\n"
00104     "                 If the transformation is specified by -matvec, then\n"
00105     "                 the matrix-vector transformation is applied to the\n"
00106     "                 (x,y,z) vectors directly, with no coordinate shifting.\n"
00107     "\n"
00108     "               * AFNI .1D files with 3 columns\n"
00109     "                   x y z\n"
00110     "                   x y z\n"
00111     "                   et cetera...\n"
00112     "                 In this case, each (x,y,z) triple is transformed and\n"
00113     "                 written to the output.  Lines that cannot be scanned\n"
00114     "                 as 3 floats are treated as comments and are passed\n"
00115     "                 through to the output unchanged.\n"
00116     "               N.B.: AFNI (x,y,z) coordinates are in DICOM order:\n"
00117     "                   -x = Right     +x = Left\n"
00118     "                   -y = Anterior  +y = Posterior\n"
00119     "                   -z = Inferior  +z = Superior\n"
00120     "\n"
00121     " -output ooo = Write the output to file 'ooo' (to stdout if 'ooo'\n"
00122     "               is '-', or if the -output option is missing).  If the\n"
00123     "               file already exists, it will not be overwritten unless\n"
00124     "               the -force option is also used.\n"
00125     "\n"
00126     " -force      = If the output file already exists, -force can be\n"
00127     "               used to overwrite it.  If you want to use -force,\n"
00128     "               it must come before -output on the command line.\n"
00129     "\n"
00130     "EXAMPLES:\n"
00131     "\n"
00132     "  Vecwarp -apar fred+tlrc -input fred.orig.coord > fred.tlrc.coord\n"
00133     "\n"
00134     "This transforms the vectors defined in original coordinates to\n"
00135     "Talairach coordinates, using the transformation previously defined\n"
00136     "by AFNI markers.\n"
00137     "\n"
00138     "  Vecwarp -apar fred+tlrc -input fred.tlrc.coord -backward > fred.test.coord\n"
00139     "\n"
00140     "This does the reverse transformation; fred.test.coord should differ from\n"
00141     "fred.orig.coord only by roundoff error.\n"
00142     "\n"
00143     "Author: RWCox - October 2001\n"
00144    ) ;
00145    exit(0) ;
00146 }