00040 {
00041    THD_3dim_dataset *inset , *outset=NULL , *newgset=NULL ;
00042    int nxin,nyin,nzin,nxyzin,nvals , ival ;
00043    int nxout,nyout,nzout,nxyzout ;
00044    char * prefix = "warped" ;
00045    int nopt=1 , verb=0 , zpad=0 , fsl=0 ;
00046    int use_matvec=0 ;
00047    float xbot,xtop , ybot,ytop , zbot,ztop ;
00048    int use_newgrid=0 ;
00049    float ddd_newgrid=0.0 , fac ;
00050    MRI_IMAGE *inim , *outim , *wim ;
00051    void *newggg=NULL ; int gflag=0 ;
00052 
00053    int tta2mni=0 , mni2tta=0 ;   /* 11 Mar 2004 */
00054    int matdefined=0 ;
00055 
00056    THD_coorder cord ;
00057    ATR_float *atr_matfor=NULL , *atr_matinv=NULL ;
00058 
00059 #if 0
00060    MRI_IMAGE *matflim=NULL ;
00061    float     *matflar=NULL ;
00062 #endif
00063 
00064    /*-- help? --*/
00065 
00066    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00067      printf("\n"
00068             "Usage: 3dWarp [options] dataset\n"
00069             "Warp (spatially transform) a 3D dataset.\n"
00070             "--------------------------\n"
00071             "Transform Defining Options: [exactly one of these must be used]\n"
00072             "--------------------------\n"
00073             "  -matvec_in2out mmm = Read a 3x4 affine transform matrix+vector\n"
00074             "                        from file 'mmm':\n"
00075             "                         x_out = Matrix x_in + Vector\n"
00076             "\n"
00077             "  -matvec_out2in mmm = Read a 3x4 affine transform matrix+vector\n"
00078             "                         from file 'mmm':\n"
00079             "                         x_in = Matrix x_out + Vector\n"
00080             "\n"
00081             "     ** N.B.: The coordinate vectors described above are\n"
00082             "               defined in DICOM ('RAI') coordinate order.\n"
00083             "               (Also see the '-fsl_matvec option, below.)\n"
00084             "     ** N.B.: Using the special name 'IDENTITY' for 'mmm'\n"
00085             "               means to use the identity matrix.\n"
00086             "     ** N.B.: You can put the matrix on the command line\n"
00087             "               directly by using an argument of the form\n"
00088             "       'MATRIX(a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34)'\n"
00089             "               in place of 'mmm', where the aij values are the\n"
00090             "               matrix entries (aij = i-th row, j-th column),\n"
00091             "               separated by commas.\n"
00092             "             * You will need the 'forward single quotes' around\n"
00093             "               the argument.\n"
00094             "\n"
00095             "  -tta2mni = Transform a dataset in Talairach-Tournoux Atlas\n"
00096             "              coordinates to MNI-152 coordinates.\n"
00097             "  -mni2tta = Transform a dataset in MNI-152 coordinates to\n"
00098             "              Talairach-Tournoux Atlas coordinates.\n"
00099             "\n"
00100             "  -matparent mset = Read in the matrix from WARPDRIVE_MATVEC_*\n"
00101             "                     attributes in the header of dataset 'mset',\n"
00102             "                     which must have been created by program\n"
00103             "                     3dWarpDrive.  In this way, you can apply\n"
00104             "                     a transformation matrix computed from\n"
00105             "                     in 3dWarpDrive to another dataset.\n"
00106             "\n"
00107             "     ** N.B.: The above option is analogous to the -rotparent\n"
00108             "                option in program 3drotate.  Use of -matparent\n"
00109             "                should be limited to datasets whose spatial\n"
00110             "                coordinate system corresponds to that which\n"
00111             "                was used for input to 3dWarpDrive (i.e., the\n"
00112             "                input to 3dWarp should overlay properly with\n"
00113             "                the input to 3dWarpDrive that generated the\n"
00114             "                -matparent dataset).\n"
00115             "              Sample usages:\n"
00116             " 3dWarpDrive -affine_general -base d1+orig -prefix d2WW -twopass -input d2+orig\n"
00117             " 3dWarp      -matparent d2WW+orig -prefix epi2WW epi2+orig\n"
00118 #if 0
00119             "\n"
00120             "  -matfile mname  = Read in the file 'mname', which consists\n"
00121             "                     of 12 ASCII-formatted numbers per line.\n"
00122             "                     The i-th input line defines a 3x4 matrix\n"
00123             "                     which is used to transform the i-th\n"
00124             "                     sub-brick of the input dataset.\n"
00125 #endif
00126             "\n"
00127             "-----------------------\n"
00128             "Other Transform Options:\n"
00129             "-----------------------\n"
00130             "  -linear     }\n"
00131             "  -cubic      } = Chooses spatial interpolation method.\n"
00132             "  -NN         } =   [default = linear]\n"
00133             "  -quintic    }\n"
00134             "\n"
00135             "  -fsl_matvec   = Indicates that the matrix file 'mmm' uses FSL\n"
00136             "                    ordered coordinates ('LPI').  For use with\n"
00137             "                    matrix files from FSL and SPM.\n"
00138             "\n"
00139             "  -newgrid ddd  = Tells program to compute new dataset on a\n"
00140             "                    new 3D grid, with spacing of 'ddd' mmm.\n"
00141             "                  * If this option is given, then the new\n"
00142             "                    3D region of space covered by the grid\n"
00143             "                    is computed by warping the 8 corners of\n"
00144             "                    the input dataset, then laying down a\n"
00145             "                    regular grid with spacing 'ddd'.\n"
00146             "                  * If this option is NOT given, then the\n"
00147             "                    new dataset is computed on the old\n"
00148             "                    dataset's grid.\n"
00149             "\n"
00150             "  -gridset ggg  = Tells program to compute new dataset on the\n"
00151             "                    same grid as dataset 'ggg'.\n"
00152             "\n"
00153             "  -zpad N       = Tells program to pad input dataset with 'N'\n"
00154             "                    planes of zeros on all sides before doing\n"
00155             "                    transformation.\n"
00156             "---------------------\n"
00157             "Miscellaneous Options:\n"
00158             "---------------------\n"
00159 #if 0
00160             "  -verb         = Print out some information along the way.\n"
00161 #endif
00162             "  -prefix ppp   = Sets the prefix of the output dataset.\n"
00163             "\n"
00164            ) ;
00165      exit(0) ;
00166    }
00167 
00168    /*-- startup mechanics --*/
00169 
00170    mainENTRY("3dWarp main"); machdep(); AFNI_logger("3dWarp",argc,argv);
00171    PRINT_VERSION("3dWarp") ;
00172 
00173    THD_coorder_fill( my_getenv("AFNI_ORIENT") , &cord ) ;  /* 12 Mar 2004 */
00174 
00175    /*-- command line options --*/
00176 
00177    while( nopt < argc && argv[nopt][0] == '-' ){
00178 
00179      /*----- 06 Jul 2005 -----*/
00180 
00181 #if 0
00182      if( strncmp(argv[nopt],"-matfile",8) == 0 ){
00183        if( matdefined )
00184          ERROR_exit("-matfile: Matrix already defined!\n");
00185        if( ++nopt >= argc )
00186          ERROR_exit("Need argument after -matfile!\n");
00187        matflim = mri_read_ascii( argv[nopt] ) ;
00188        if( matflim == NULL )
00189          ERROR_exit("Can't read from -matfile %s\n",argv[nopt]);
00190        if( matflim->nx < 12 )
00191          ERROR_exit("Need 12 numbers per line in -matfile %s\n",argv[nopt]) ;
00192 
00193        matflar = MRI_FLOAT_PTR(matflim) ;
00194        matdefined = 1 ; nopt++ ; continue ;
00195      }
00196 #endif
00197 
00198      /*-----*/
00199 
00200      if( strncmp(argv[nopt],"-matparent",10) == 0 ){
00201        ATR_float *atr ; float *matar ; THD_3dim_dataset *matparset=NULL ;
00202 
00203        if( matdefined )
00204          ERROR_exit("-matparent: Matrix already defined!\n");
00205        if( ++nopt >= argc )
00206          ERROR_exit("Need argument after -matparent!\n");
00207        matparset = THD_open_dataset( argv[nopt] ) ;
00208        if( matparset == NULL )
00209          ERROR_exit("Can't open -matparent %s\n",argv[nopt]);
00210        atr = THD_find_float_atr( matparset->dblk, "WARPDRIVE_MATVEC_INV_000000" );
00211        if( atr != NULL ) atr_matinv = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
00212        atr = THD_find_float_atr( matparset->dblk, "WARPDRIVE_MATVEC_FOR_000000" );
00213        if( atr != NULL ) atr_matfor = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
00214        if( atr_matinv == NULL ||  atr_matinv->nfl < 12 )
00215          ERROR_exit("-matparent %s doesn't have WARPDRIVE attributes!?\n",argv[nopt]) ;
00216 
00217        matar = atr_matinv->fl ;
00218        if( strstr(argv[nopt-1],"INV") == NULL ){
00219          LOAD_MAT  ( dicom_in2out.mm, matar[0],matar[1],matar[2],
00220                                       matar[4],matar[5],matar[6],
00221                                       matar[8],matar[9],matar[10] ) ;
00222          LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;
00223          dicom_out2in = INV_VECMAT( dicom_in2out ) ;
00224        } else {
00225          LOAD_MAT  ( dicom_out2in.mm, matar[0],matar[1],matar[2],
00226                                       matar[4],matar[5],matar[6],
00227                                       matar[8],matar[9],matar[10] ) ;
00228          LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;
00229          dicom_in2out = INV_VECMAT( dicom_out2in ) ;
00230        }
00231 
00232        DSET_delete(matparset) ;
00233 
00234        use_matvec = MATVEC_FOR ; matdefined = 1 ; nopt++ ; continue ;
00235      }
00236 
00237      /*----- 11 Mar 2004 -----*/
00238 
00239      if( strcmp(argv[nopt],"-tta2mni") == 0 ){
00240        if( matdefined )
00241          ERROR_exit("-tta2mni: Matrix already defined!\n");
00242        matdefined = 1 ; tta2mni = 1 ; nopt++ ; continue ;
00243      }
00244 
00245      if( strcmp(argv[nopt],"-mni2tta") == 0 ){
00246        if( matdefined )
00247          ERROR_exit("-mni2tta: Matrix already defined!\n");
00248        matdefined = 1 ; mni2tta = 1 ; nopt++ ; continue ;
00249      }
00250 
00251      /*-----*/
00252 
00253      if( strcmp(argv[nopt],"-fsl_matvec") == 0 ){
00254        fsl = 1 ; nopt++ ; continue ;
00255      }
00256 
00257      /*-----*/
00258 
00259      if( strcmp(argv[nopt],"-gridset") == 0 ){
00260        if( use_newgrid )
00261          ERROR_exit("Can't use -gridset with -newgrid!\n");
00262        if( newgset != NULL )
00263          ERROR_exit("Can't use -gridset twice!\n");
00264        if( ++nopt >= argc )
00265          ERROR_exit("Need argument after -gridset!\n");
00266        newgset = THD_open_dataset( argv[nopt] ) ;
00267        if( newgset == NULL )
00268          ERROR_exit("Can't open -gridset %s\n",argv[nopt]);
00269        nopt++ ; continue ;
00270      }
00271 
00272      /*-----*/
00273 
00274      if( strcmp(argv[nopt],"-zpad") == 0 ){
00275        if( ++nopt >= argc )
00276          ERROR_exit("Need argument after -zpad!\n");
00277        zpad = (int) strtod( argv[nopt] , NULL ) ;
00278        if( zpad < 0 )
00279          ERROR_exit("Illegal negative argument after -zpad!\n");
00280        else if( zpad == 0 )
00281          INFO_message("NOTA BENE: -zpad is 0 ==> ignored\n") ;
00282        nopt++ ; continue ;
00283      }
00284 
00285      /*-----*/
00286 
00287      if( strcmp(argv[nopt],"-newgrid") == 0 ){
00288        if( newgset != NULL )
00289          ERROR_exit("Can't use -newgrid with -gridset!\n");
00290        if( use_newgrid )
00291          ERROR_exit("Can't use -newgrid twice!\n");
00292        if( ++nopt >= argc )
00293          ERROR_exit("Need argument after -newgrid!\n");
00294 
00295        ddd_newgrid = strtod( argv[nopt] , NULL ) ;
00296        if( ddd_newgrid <= 0.0 )
00297          ERROR_exit("Illegal argument after -newgrid!\n");
00298 
00299        use_newgrid = 1 ; nopt++ ; continue ;
00300      }
00301 
00302      /*-----*/
00303 
00304      if( strcmp(argv[nopt],"-NN")     == 0 ){
00305        mri_warp3D_method( MRI_NN )     ; nopt++ ; continue ;
00306      }
00307      if( strcmp(argv[nopt],"-linear") == 0 ){
00308        mri_warp3D_method( MRI_LINEAR ) ; nopt++ ; continue ;
00309      }
00310      if( strcmp(argv[nopt],"-cubic")  == 0 ){
00311        mri_warp3D_method( MRI_CUBIC )  ; nopt++ ; continue ;
00312      }
00313      if( strcmp(argv[nopt],"-quintic") == 0 ){
00314        mri_warp3D_method( MRI_QUINTIC ); nopt++ ; continue ;  /* 06 Aug 2003 */
00315      }
00316 
00317      /*-----*/
00318 
00319      if( strcmp(argv[nopt],"-prefix") == 0 ){
00320        if( ++nopt >= argc )
00321          ERROR_exit("Need an argument after -prefix!\n");
00322        if( !THD_filename_ok(argv[nopt]) )
00323          ERROR_exit("-prefix argument is invalid!\n");
00324 
00325        prefix = argv[nopt] ; nopt++ ; continue ;
00326      }
00327 
00328      /*-----*/
00329 
00330      if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00331        verb++ ; nopt++ ; continue ;
00332      }
00333 
00334      /*-----*/
00335 
00336      if( strncmp(argv[nopt],"-matvec_",8) == 0 ){
00337        MRI_IMAGE *matim ; float *matar , dm ;
00338 
00339        if( matdefined )
00340          ERROR_exit("%s: matrix already defined!\n",argv[nopt]);
00341        if( ++nopt >= argc )
00342          ERROR_exit("Need an argument after -matvec!\n");
00343        if( strcmp(argv[nopt],"IDENTITY") == 0 ){ /* load identity matrix */
00344          matim = mri_new( 4,3 , MRI_float ) ;    /* will be all zero */
00345          matar = MRI_FLOAT_PTR(matim) ;
00346          matar[0] = matar[5] = matar[10] = 1.0 ;
00347        } else if( strncmp(argv[nopt],"MATRIX(",7) == 0 ){ /* matrix from arg */
00348          int nn ;
00349          matim = mri_new( 4,3 , MRI_float ) ;    /* will be all zero */
00350          matar = MRI_FLOAT_PTR(matim) ;
00351          nn = sscanf(argv[nopt],"MATRIX(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f",
00352                      matar+0 , matar+1 , matar+2 , matar+3 ,
00353                      matar+4 , matar+5 , matar+6 , matar+7 ,
00354                      matar+8 , matar+9 , matar+10, matar+11 ) ;
00355          if( nn < 12 )
00356            ERROR_exit("Could only decode %d numbers from %s\n",nn,argv[nopt]) ;
00357        } else {
00358          matim = mri_read_ascii( argv[nopt] ) ;
00359          if( matim == NULL )
00360            ERROR_exit("Can't read -matvec file %s\n",argv[nopt]);
00361          if( matim->nx != 4 || matim->ny < 3 )
00362            ERROR_exit("-matvec file not a 3x4 matrix!\n");
00363          matar = MRI_FLOAT_PTR(matim) ;
00364        }
00365 
00366        use_matvec = (strstr(argv[nopt-1],"_in2out") != NULL) ? MATVEC_FOR : MATVEC_BAC ;
00367 
00368        switch( use_matvec ){
00369 
00370          case MATVEC_FOR:
00371            LOAD_MAT  ( dicom_in2out.mm, matar[0],matar[1],matar[2],
00372                                         matar[4],matar[5],matar[6],
00373                                         matar[8],matar[9],matar[10] ) ;
00374            LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;
00375 
00376            dm = MAT_DET( dicom_in2out.mm ) ;
00377            if( dm == 0.0 )
00378              ERROR_exit("Determinant of matrix is 0\n");
00379 
00380            dicom_out2in = INV_VECMAT( dicom_in2out ) ;
00381          break ;
00382 
00383          case MATVEC_BAC:
00384            LOAD_MAT  ( dicom_out2in.mm, matar[0],matar[1],matar[2],
00385                                         matar[4],matar[5],matar[6],
00386                                         matar[8],matar[9],matar[10] ) ;
00387            LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;
00388 
00389            dm = MAT_DET( dicom_out2in.mm ) ;
00390            if( dm == 0.0 )
00391              ERROR_exit("Determinant of matrix is 0\n");
00392 
00393            dicom_in2out = INV_VECMAT( dicom_out2in ) ;
00394          break ;
00395        }
00396 
00397        matdefined = 1 ; nopt++ ; continue ;
00398      }
00399 
00400      /*-----*/
00401 
00402      ERROR_exit("Unknown option %s\n",argv[nopt]) ;
00403 
00404    } /* end of loop over options */
00405 
00406    /*-- check to see if we have a warp specified --*/
00407 
00408    if( !matdefined )
00409      ERROR_exit("No transformation on command line!?\n");
00410 
00411    if( tta2mni || mni2tta ) fsl = 0 ;
00412 
00413    /*-- 06 Aug 2003: if FSL matrix, flip stuff around --*/
00414 
00415    if( use_matvec && fsl ){
00416      dicom_in2out.mm.mat[0][2] = -dicom_in2out.mm.mat[0][2] ;
00417      dicom_in2out.mm.mat[1][2] = -dicom_in2out.mm.mat[1][2] ;
00418      dicom_in2out.mm.mat[2][0] = -dicom_in2out.mm.mat[2][0] ;
00419      dicom_in2out.mm.mat[2][1] = -dicom_in2out.mm.mat[2][1] ;
00420      dicom_in2out.vv.xyz[0]    = -dicom_in2out.vv.xyz[0]    ;
00421      dicom_in2out.vv.xyz[1]    = -dicom_in2out.vv.xyz[1]    ;
00422 
00423      dicom_out2in.mm.mat[0][2] = -dicom_out2in.mm.mat[0][2] ;
00424      dicom_out2in.mm.mat[1][2] = -dicom_out2in.mm.mat[1][2] ;
00425      dicom_out2in.mm.mat[2][0] = -dicom_out2in.mm.mat[2][0] ;
00426      dicom_out2in.mm.mat[2][1] = -dicom_out2in.mm.mat[2][1] ;
00427      dicom_out2in.vv.xyz[0]    = -dicom_out2in.vv.xyz[0]    ;
00428      dicom_out2in.vv.xyz[1]    = -dicom_out2in.vv.xyz[1]    ;
00429    }
00430 
00431    /*-- 1 remaining argument should be a dataset --*/
00432 
00433    if( nopt != argc-1 ){
00434      ERROR_message("Command line should have exactly 1 dataset!") ;
00435      ERROR_exit   ("Whereas there seems to be %d of them!\n", argc-nopt ) ;
00436    }
00437 
00438    /*-- input dataset header --*/
00439 
00440    inset = THD_open_dataset( argv[nopt] ) ;
00441    if( !ISVALID_DSET(inset) )
00442      ERROR_exit("Can't open dataset %s\n",argv[nopt]);
00443 
00444    /*- 11 Mar 2004: check for coherency -*/
00445 
00446    if( tta2mni || mni2tta ){
00447      if( inset->view_type != VIEW_TALAIRACH_TYPE ){
00448        WARNING_message("Input dataset not in +tlrc view!\n") ;
00449      } else {
00450        ATR_string *ast ;
00451        ast = THD_find_string_atr( inset->dblk , "TLRC_SPACE" ) ;
00452        if( ast != NULL ){
00453          if( strcmp(ast->ch,"MNI-152") == 0 && tta2mni ){
00454            WARNING_message("Input dataset appears to be MNI-152 already!\n");
00455          }
00456        }
00457      }
00458    }
00459 
00460    /*-- do the heavy lifting --*/
00461 
00462    if( use_newgrid ){
00463      newggg = &ddd_newgrid ; gflag = WARP3D_NEWGRID ;
00464    } else if( newgset != NULL ){
00465      newggg = newgset      ; gflag = WARP3D_NEWDSET ;
00466    }
00467 
00468    if( use_matvec ){
00469      outset = THD_warp3D_affine( inset , dicom_out2in ,
00470                                  newggg , prefix , zpad , gflag ) ;
00471    } else if( tta2mni ){
00472      outset = THD_warp3D_tta2mni( inset , newggg , prefix,zpad,gflag ) ;
00473    } else if( mni2tta ){
00474      outset = THD_warp3D_mni2tta( inset , newggg , prefix,zpad,gflag ) ;
00475    }
00476 
00477    if( outset == NULL )
00478      ERROR_exit("THD_warp3D() fails for some reason!\n") ;
00479 
00480    /*-- polish up the new dataset info --*/
00481 
00482    tross_Copy_History( inset , outset ) ;
00483    tross_Make_History( "3dWarp" , argc,argv , outset ) ;
00484 
00485    /*- 11 Mar 2004: mark output dataset as being in MNI space -*/
00486 
00487    if( tta2mni ){
00488      THD_set_string_atr( outset->dblk , "TLRC_SPACE" , "MNI-152" ) ;
00489    }
00490 
00491    /*- 03 Aug 2005: save WARPDRIVE attributes, if present -*/
00492 
00493    if( atr_matfor != NULL ) THD_insert_atr( outset->dblk, (ATR_any *)atr_matfor );
00494    if( atr_matinv != NULL ) THD_insert_atr( outset->dblk, (ATR_any *)atr_matinv );
00495 
00496    /*- actually do the writositing -*/
00497 
00498    DSET_delete( inset ) ;
00499    DSET_write( outset ) ; if( verb ) WROTE_DSET(outset) ;
00500    exit(0) ;
00501 }

00019 {
00020    THD_fvec3 xxx ;
00021    LOAD_FVEC3( xxx , xin,yin,zin ) ;
00022    xxx = VECMAT_VEC( dicom_in2out , xxx ) ;
00023    UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
00024 }

00030 {
00031    THD_fvec3 xxx ;
00032    LOAD_FVEC3( xxx , xin,yin,zin ) ;
00033    xxx = VECMAT_VEC( dicom_out2in , xxx ) ;
00034    UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
00035 }