00001
00002
00003
00004
00005
00006
00007 #include "thd_shear3d.h"
00008
00009
00010
00011 static void rotate_stdin_points( THD_dfvec3, THD_dmat33, int,THD_dfvec3 ) ;
00012
00013
00014
00015 #define MATVEC_DICOM 1
00016 #define MATVEC_ORDER 2
00017
00018 static int verb=0 ;
00019
00020
00021
00022 #define MODE_DFILE 1
00023 #define MODE_1DFILE 2
00024 MRI_IMAGE * get_dfile_params( char * fname , int mode ) ;
00025
00026
00027
00028 int main( int argc , char * argv[] )
00029 {
00030 THD_3dim_dataset * dset=NULL ;
00031 char * new_prefix = "rota" , * cpt ;
00032 float dx=0 , dy=0 , dz=0 ;
00033 int ax1=0,ax2=1,ax3=2 , adx,ady,adz ;
00034 char cdx,cdy,cdz ;
00035 float th1=0.0,th2=0.0,th3=0.0 ;
00036 int iopt , nvox , rotarg=-1 , dcode=-1 , ival,nval , ihand ;
00037 float * fvol ;
00038 double cputim ;
00039 int clipit=1 ;
00040 float cbot,ctop ;
00041
00042 int matvec=0 ;
00043 THD_dmat33 rmat , pp,ppt ;
00044 THD_dfvec3 tvec ;
00045
00046 int dopoints=0 , doorigin=0 ;
00047 double xo=0.0 , yo=0.0 , zo=0.0 ;
00048
00049 THD_3dim_dataset *rotpar_dset=NULL , *gridpar_dset=NULL ;
00050
00051 int zpad=0 ;
00052
00053 char *dname=NULL ;
00054 int dmode=0 ;
00055 MRI_IMAGE *dim=NULL ;
00056 float *dar=NULL ;
00057 int ndar=0 ;
00058 int skipit=0 ;
00059
00060
00061
00062 LOAD_DIAG_DMAT(rmat,1.0,1.0,1.0) ;
00063
00064
00065
00066 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){
00067 printf(
00068 "Usage: 3drotate [options] dataset\n"
00069 "Rotates and/or translates all bricks from an AFNI dataset.\n"
00070 "'dataset' may contain a sub-brick selector list.\n"
00071 "\n"
00072 "GENERIC OPTIONS:\n"
00073 " -prefix fname = Sets the output dataset prefix name to be 'fname'\n"
00074 " -verbose = Prints out progress reports (to stderr)\n"
00075 "\n"
00076 "OPTIONS TO SPECIFY THE ROTATION/TRANSLATION:\n"
00077 "-------------------------------------------\n"
00078 "*** METHOD 1 = direct specification:\n"
00079 "At most one of these shift options can be used:\n"
00080 " -ashift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n"
00081 " AFTER rotation.\n"
00082 " -bshift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n"
00083 " BEFORE rotation.\n"
00084 " The shift distances by default are along the (x,y,z) axes of the dataset\n"
00085 " storage directions (see the output of '3dinfo dataset'). To specify them\n"
00086 " anatomically, you can suffix a distance with one of the symbols\n"
00087 " 'R', 'L', 'A', 'P', 'I', and 'S', meaning 'Right', 'Left', 'Anterior',\n"
00088 " 'Posterior', 'Inferior', and 'Superior', respectively.\n"
00089 "\n"
00090 " -rotate th1 th2 th3\n"
00091 " Specifies the 3D rotation to be composed of 3 planar rotations:\n"
00092 " 1) 'th1' degrees about the 1st axis, followed by\n"
00093 " 2) 'th2' degrees about the (rotated) 2nd axis, followed by\n"
00094 " 3) 'th3' degrees about the (doubly rotated) 3rd axis.\n"
00095 " Which axes are used for these rotations is specified by placing\n"
00096 " one of the symbols 'R', 'L', 'A', 'P', 'I', and 'S' at the end\n"
00097 " of each angle (e.g., '10.7A'). These symbols denote rotation\n"
00098 " about the 'Right-to-Left', 'Left-to-Right', 'Anterior-to-Posterior',\n"
00099 " 'Posterior-to-Anterior', 'Inferior-to-Superior', and\n"
00100 " 'Superior-to-Inferior' axes, respectively. A positive rotation is\n"
00101 " defined by the right-hand rule.\n"
00102 "\n"
00103 "*** METHOD 2 = copy from output of 3dvolreg:\n"
00104 " -rotparent rset\n"
00105 " Specifies that the rotation and translation should be taken from the\n"
00106 " first 3dvolreg transformation found in the header of dataset 'rset'.\n"
00107 " -gridparent gset\n"
00108 " Specifies that the output dataset of 3drotate should be shifted to\n"
00109 " match the grid of dataset 'gset'. Can only be used with -rotparent.\n"
00110 " This dataset should be one this is properly aligned with 'rset' when\n"
00111 " overlaid in AFNI.\n"
00112 " * If -rotparent is used, then don't use -matvec, -rotate, or -[ab]shift.\n"
00113 " * If 'gset' has a different number of slices than the input dataset,\n"
00114 " then the output dataset will be zero-padded in the slice direction\n"
00115 " to match 'gset'.\n"
00116 " * These options are intended to be used to align datasets between sessions:\n"
00117 " S1 = SPGR from session 1 E1 = EPI from session 1\n"
00118 " S2 = SPGR from session 2 E2 = EPI from session 2\n"
00119 " 3dvolreg -twopass -twodup -base S1+orig -prefix S2reg S2+orig\n"
00120 " 3drotate -rotparent S2reg+orig -gridparent E1+orig -prefix E2reg E2+orig\n"
00121 " The result will have E2reg rotated from E2 in the same way that S2reg\n"
00122 " was from S2, and also shifted/padded (as needed) to overlap with E1.\n"
00123 "\n"
00124 "*** METHOD 3 = give the transformation matrix/vector directly:\n"
00125 " -matvec_dicom mfile\n"
00126 " -matvec_order mfile\n"
00127 " Specifies that the rotation and translation should be read from file\n"
00128 " 'mfile', which should be in the format\n"
00129 " u11 u12 u13 v1\n"
00130 " u21 u22 u23 v2\n"
00131 " u31 u32 u33 u3\n"
00132 " where each 'uij' and 'vi' is a number. The 3x3 matrix [uij] is the\n"
00133 " orthogonal matrix of the rotation, and the 3-vector [vi] is the -ashift\n"
00134 " vector of the translation.\n"
00135 "\n"
00136 "*** METHOD 4 = copy the transformation from 3dTagalign:\n"
00137 " -matvec_dset mset\n"
00138 " Specifies that the rotation and translation should be read from\n"
00139 " the .HEAD file of dataset 'mset', which was created by program\n"
00140 " 3dTagalign.\n"
00141 " * If -matvec_dicom is used, the matrix and vector are given in Dicom\n"
00142 " coordinate order (+x=L, +y=P, +z=S). This is the option to use\n"
00143 " if mfile is generated using 3dTagalign -matvec mfile.\n"
00144 " * If -matvec_order is used, the the matrix and vector are given in the\n"
00145 " coordinate order of the dataset axes, whatever they may be.\n"
00146 " * You can't mix -matvec_* options with -rotate and -*shift.\n"
00147 "\n"
00148 "*** METHOD 5 = input rotation+shift parameters from an ASCII file:\n"
00149 " -dfile dname *OR* -1Dfile dname\n"
00150 " With these methods, the movement parameters for each sub-brick\n"
00151 " of the input dataset are read from the file 'dname'. This file\n"
00152 " should consist of columns of numbers in ASCII format. Six (6)\n"
00153 " numbers are read from each line of the input file. If the\n"
00154 " '-dfile' option is used, each line of the input should be at\n"
00155 " least 7 numbers, and be of the form\n"
00156 " ignored roll pitch yaw dS dL dP\n"
00157 " If the '-1Dfile' option is used, then each line of the input\n"
00158 " should be at least 6 numbers, and be of the form\n"
00159 " roll pitch yaw dS dL dP\n"
00160 " (These are the forms output by the '-dfile' and\n"
00161 " '-1Dfile' options of program 3dvolreg; see that\n"
00162 " program's -help output for the hideous details.)\n"
00163 " The n-th sub-brick of the input dataset will be transformed\n"
00164 " using the parameters from the n-th line of the dname file.\n"
00165 " If the dname file doesn't contain as many lines as the\n"
00166 " input dataset has sub-bricks, then the last dname line will\n"
00167 " be used for all subsequent sub-bricks. Excess columns or\n"
00168 " rows will be ignored.\n"
00169 " N.B.: Rotation is always about the center of the volume.\n"
00170 " If the parameters are derived from a 3dvolreg run\n"
00171 " on a dataset with a different center in xyz-space,\n"
00172 " the results may not be what you want!\n"
00173 " N.B.: You can't use -dfile/-1Dfile with -points (infra).\n"
00174 "\n"
00175 "POINTS OPTIONS (instead of datasets):\n"
00176 "------------------------------------\n"
00177 " -points\n"
00178 " -origin xo yo zo\n"
00179 " These options specify that instead of rotating a dataset, you will\n"
00180 " be rotating a set of (x,y,z) points. The points are read from stdin.\n"
00181 " * If -origin is given, the point (xo,yo,zo) is used as the center for\n"
00182 " the rotation.\n"
00183 " * If -origin is NOT given, and a dataset is given at the end of the\n"
00184 " command line, then the center of the dataset brick is used as\n"
00185 " (xo,yo,zo). The dataset will NOT be rotated if -points is given.\n"
00186 " * If -origin is NOT given, and NO dataset is given at the end of the\n"
00187 " command line, then xo=yo=zo=0 is assumed. You probably don't\n"
00188 " want this.\n"
00189 " * (x,y,z) points are read from stdin as 3 ASCII-formatted numbers per\n"
00190 " line, as in 3dUndump. Any succeeding numbers on input lines will\n"
00191 " be copied to the output, which will be written to stdout.\n"
00192 " * The input (x,y,z) coordinates are taken in the same order as the\n"
00193 " axes of the input dataset. If there is no input dataset, then\n"
00194 " negative x = R positive x = L }\n"
00195 " negative y = A positive y = P } e.g., the DICOM order\n"
00196 " negative z = I positive z = S }\n"
00197 " One way to dump some (x,y,z) coordinates from a dataset is:\n"
00198 "\n"
00199 " 3dmaskdump -mask something+tlrc -o xyzfilename -noijk\n"
00200 " '3dcalc( -a dset+tlrc -expr x -datum float )'\n"
00201 " '3dcalc( -a dset+tlrc -expr y -datum float )'\n"
00202 " '3dcalc( -a dset+tlrc -expr z -datum float )'\n"
00203 "\n"
00204 " (All of this should be on one command line.)\n"
00205 "============================================================================\n"
00206 "\n"
00207 "Example: 3drotate -prefix Elvis -bshift 10S 0 0 -rotate 30R 0 0 Sinatra+orig\n"
00208 "\n"
00209 "This will shift the input 10 mm in the superior direction, followed by a 30\n"
00210 "degree rotation about the Right-to-Left axis (i.e., nod the head forward).\n"
00211 "\n"
00212 "============================================================================\n"
00213 "Algorithm: The rotation+shift is decomposed into 4 1D shearing operations\n"
00214 " (a 3D generalization of Paeth's algorithm). The interpolation\n"
00215 " (i.e., resampling) method used for these shears can be controlled\n"
00216 " by the following options:\n"
00217 "\n"
00218 " -Fourier = Use a Fourier method (the default: most accurate; slowest).\n"
00219 " -NN = Use the nearest neighbor method.\n"
00220 " -linear = Use linear (1st order polynomial) interpolation (least accurate).\n"
00221 " -cubic = Use the cubic (3rd order) Lagrange polynomial method.\n"
00222 " -quintic = Use the quintic (5th order) Lagrange polynomial method.\n"
00223 " -heptic = Use the heptic (7th order) Lagrange polynomial method.\n"
00224 "\n"
00225 " -Fourier_nopad = Use the Fourier method WITHOUT padding\n"
00226 " * If you don't mind - or even want - the wraparound effect\n"
00227 " * Works best if dataset grid size is a power of 2, possibly\n"
00228 " times powers of 3 and 5, in all directions being altered.\n"
00229 " * The main use would seem to be to un-wraparound poorly\n"
00230 " reconstructed images, by using a shift; for example:\n"
00231 " 3drotate -ashift 30A 0 0 -Fourier_nopad -prefix Anew A+orig\n"
00232 " * This option is also available in the Nudge Dataset plugin.\n"
00233 "\n"
00234 " -clipit = Clip results to input brick range [now the default].\n"
00235 " -noclip = Don't clip results to input brick range.\n"
00236 "\n"
00237 " -zpad n = Zeropad around the edges by 'n' voxels during rotations\n"
00238 " (these edge values will be stripped off in the output)\n"
00239 " N.B.: Unlike to3d, in this program '-zpad' adds zeros in\n"
00240 " all directions.\n"
00241 " N.B.: The environment variable AFNI_ROTA_ZPAD can be used\n"
00242 " to set a nonzero default value for this parameter.\n"
00243 ) ;
00244
00245 printf("\n" MASTER_SHORTHELP_STRING ) ;
00246 exit(0) ;
00247 }
00248
00249 mainENTRY("3drotate main") ; machdep() ; PRINT_VERSION("3drotate") ;
00250
00251
00252
00253 { int new_argc ; char ** new_argv ;
00254 addto_args( argc , argv , &new_argc , &new_argv ) ;
00255 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00256 }
00257
00258 AFNI_logger("3drotate",argc,argv) ;
00259
00260 #define ERREX(str) (fprintf(stderr,"*** %s\n",str),exit(1))
00261
00262 iopt = 1 ;
00263 while( iopt < argc && argv[iopt][0] == '-' ){
00264
00265 if( strncmp(argv[iopt],"-zpad",5) == 0 ){
00266 if( zpad > 0 )
00267 fprintf(stderr,"+++ WARNING: second -zpad option!\n") ;
00268 zpad = (int) strtod( argv[++iopt] , NULL ) ;
00269 if( zpad < 0 ){
00270 fprintf(stderr,"*** ERROR: Can't use -zpad %d\n",zpad) ;
00271 exit(1) ;
00272 }
00273 THD_rota_setpad(zpad,zpad,zpad) ;
00274 iopt++ ; continue ;
00275 }
00276
00277 if( strncmp(argv[iopt],"-points",7) == 0 ){
00278 dopoints = 1 ;
00279 iopt++ ; continue ;
00280 }
00281
00282 if( strncmp(argv[iopt],"-origin",7) == 0 ){
00283 xo = strtod( argv[++iopt] , NULL ) ;
00284 yo = strtod( argv[++iopt] , NULL ) ;
00285 zo = strtod( argv[++iopt] , NULL ) ;
00286 doorigin = 1 ;
00287 iopt++ ; continue ;
00288 }
00289
00290 if( strncmp(argv[iopt],"-rotpar",7) == 0 ){
00291
00292 ATR_float *atr ;
00293
00294 if( rotpar_dset != NULL )
00295 ERREX("*** Can't use 2 -rotparent options!") ;
00296 if( matvec )
00297 ERREX("*** Can't combine -rotparent with -matvec!") ;
00298 if( dcode > 0 || rotarg > 0 )
00299 ERREX("*** Can't use -rotparent with -shift or -rotate options!") ;
00300 if( dname != NULL )
00301 ERREX("*** Can't use -rotparent with -dfile or -1Dfile options!") ;
00302
00303 rotpar_dset = THD_open_one_dataset( argv[++iopt] ) ;
00304 if( rotpar_dset == NULL )
00305 ERREX("*** Can't open -rotparent dataset!\n") ;
00306
00307 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ;
00308 if( atr == NULL || atr->nfl < 12 )
00309 ERREX("*** -rotparent dataset doesn't have VOLREG attributes!?") ;
00310
00311 iopt++ ; continue ;
00312 }
00313
00314 if( strncmp(argv[iopt],"-gridpar",7) == 0 ){
00315
00316 if( gridpar_dset != NULL )
00317 ERREX("*** Can't use -2 -gridparent options!") ;
00318
00319 gridpar_dset = THD_open_one_dataset( argv[++iopt] ) ;
00320 if( gridpar_dset == NULL )
00321 ERREX("*** Can't open -gridparent dataset!\n") ;
00322
00323 iopt++ ; continue ;
00324 }
00325
00326 if( strncmp(argv[iopt],"-matvec_",8) == 0 ){
00327
00328 MRI_IMAGE *matim ; float *matar , sum ;
00329
00330 if( matvec )
00331 ERREX("*** Can't use 2 -matvec options!") ;
00332 if( dcode > 0 || rotarg > 0 )
00333 ERREX("*** Can't use -matvec with -shift or -rotate options!") ;
00334 if( rotpar_dset != NULL )
00335 ERREX("*** Can't use -matvec with -rotparent option!") ;
00336 if( dname != NULL )
00337 ERREX("*** Can't use -matvec with -dfile or -1Dfile options!") ;
00338
00339 if( strcmp(argv[iopt],"-matvec_order") == 0 ) matvec = MATVEC_ORDER ;
00340 else matvec = MATVEC_DICOM ;
00341
00342 if( strcmp(argv[iopt],"-matvec_dset") == 0){
00343 THD_3dim_dataset * mvset ; ATR_float * atr ;
00344
00345 mvset = THD_open_dataset( argv[++iopt] ) ;
00346 if( mvset == NULL ) ERREX("*** Can't read -matvec_dset dataset!") ;
00347 atr = THD_find_float_atr( mvset->dblk , "TAGALIGN_MATVEC" ) ;
00348 if( atr == NULL || atr->nfl < 12 )
00349 ERREX("*** -matvec_dset doesn't have matrix+vector in .HEAD!") ;
00350 matar = atr->fl ;
00351 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
00352 matar[4],matar[5],matar[6],
00353 matar[8],matar[9],matar[10] ) ;
00354 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
00355 DSET_delete(mvset) ;
00356 } else {
00357 matim = mri_read_ascii( argv[++iopt] ) ;
00358 if( matim == NULL ) ERREX("Can't read -matvec file!") ;
00359 if( matim->nx != 4 || matim->ny != 3 ) ERREX("-matvec file not 4x3!") ;
00360
00361 matar = MRI_FLOAT_PTR(matim) ;
00362 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
00363 matar[4],matar[5],matar[6],
00364 matar[8],matar[9],matar[10] ) ;
00365 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
00366
00367 mri_free(matim) ;
00368 }
00369
00370
00371
00372
00373 pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ;
00374 sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0])
00375 +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1])
00376 +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0);
00377 if( sum > 0.01 ) ERREX("-matvec matrix not orthogonal!") ;
00378
00379 iopt++ ; continue ;
00380 }
00381
00382 if( strncmp(argv[iopt],"-clipit",4) == 0 ){
00383 fprintf(stderr,"++ Notice: -clipit is now the default\n") ;
00384 clipit = 1 ;
00385 iopt++ ; continue ;
00386 }
00387
00388 if( strncmp(argv[iopt],"-noclip",4) == 0 ){
00389 clipit = 0 ;
00390 iopt++ ; continue ;
00391 }
00392
00393 if( strncmp(argv[iopt],"-prefix",4) == 0 ){
00394 new_prefix = argv[++iopt] ;
00395 iopt++ ; continue ;
00396 }
00397
00398 if( strncmp(argv[iopt],"-verbose",5) == 0 ){
00399 verb = 1 ;
00400 iopt++ ; continue ;
00401 }
00402
00403 if( strcmp(argv[iopt],"-Fourier_nopad") == 0 ){
00404 THD_rota_method( MRI_FOURIER_NOPAD ) ;
00405 iopt++ ; continue ;
00406 }
00407 if( strncmp(argv[iopt],"-Fourier",4) == 0 || strncmp(argv[iopt],"-fourier",4) == 0 ){
00408 THD_rota_method( MRI_FOURIER ) ;
00409 iopt++ ; continue ;
00410 }
00411
00412 if( strncmp(argv[iopt],"-cubic",4) == 0 || strncmp(argv[iopt],"-Cubic",4) == 0 ){
00413 THD_rota_method( MRI_CUBIC ) ;
00414 iopt++ ; continue ;
00415 }
00416
00417 if( strncmp(argv[iopt],"-quintic",4) == 0 || strncmp(argv[iopt],"-Quintic",4) == 0 ){
00418 THD_rota_method( MRI_QUINTIC ) ;
00419 iopt++ ; continue ;
00420 }
00421
00422 if( strncmp(argv[iopt],"-heptic",4) == 0 || strncmp(argv[iopt],"-Heptic",4) == 0 ){
00423 THD_rota_method( MRI_HEPTIC ) ;
00424 iopt++ ; continue ;
00425 }
00426
00427 if( strncmp(argv[iopt],"-linear",4) == 0 || strncmp(argv[iopt],"-Linear",4) == 0 ){
00428 THD_rota_method( MRI_LINEAR ) ; clipit = 0 ;
00429 iopt++ ; continue ;
00430 }
00431
00432 if( strncmp(argv[iopt],"-nn",3) == 0 || strncmp(argv[iopt],"-NN",4) == 0 ){
00433 THD_rota_method( MRI_NN ) ; clipit = 0 ;
00434 iopt++ ; continue ;
00435 }
00436
00437 if( strncmp(argv[iopt],"-ashift",4) == 0 ){
00438 if( matvec ) ERREX("*** Can't use -ashift with -matvec!") ;
00439 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
00440 if( rotpar_dset != NULL )
00441 ERREX("*** Can't use -ashift with -rotparent!") ;
00442 if( dname != NULL )
00443 ERREX("*** Can't use -ashift with -dfile or -1Dfile options!") ;
00444 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
00445 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
00446 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
00447 dcode = DELTA_AFTER ;
00448 iopt++ ; continue ;
00449 }
00450
00451 if( strncmp(argv[iopt],"-bshift",4) == 0 ){
00452 if( matvec ) ERREX("*** Can't use -bshift with -matvec!") ;
00453 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
00454 if( rotpar_dset != NULL )
00455 ERREX("*** Can't use -bshift with -rotparent!") ;
00456 if( dname != NULL )
00457 ERREX("*** Can't use -bshift with -dfile or -1Dfile options!") ;
00458 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
00459 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
00460 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
00461 dcode = DELTA_BEFORE ;
00462 iopt++ ; continue ;
00463 }
00464
00465 #if 0
00466 if( strncmp(argv[iopt],"-cshift",4) == 0 ){
00467 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
00468 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
00469 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
00470 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
00471 dcode = DELTA_FIXED ;
00472 iopt++ ; continue ;
00473 }
00474 #endif
00475
00476 if( strncmp(argv[iopt],"-rotate",4) == 0 ){
00477 if( matvec ) ERREX("*** Can't use -rotate with -matvec!") ;
00478 if( rotarg > 0 ) ERREX("*** Can't have 2 -rotate options!\n") ;
00479 if( rotpar_dset != NULL )
00480 ERREX("*** Can't use -rotate with -rotparent!") ;
00481 if( dname != NULL )
00482 ERREX("*** Can't use -rotate with -dfile or -1Dfile options!") ;
00483 rotarg = iopt ;
00484 iopt += 4 ; continue ;
00485 }
00486
00487
00488
00489 if( strcmp(argv[iopt],"-dfile")==0 || strcmp(argv[iopt],"-1Dfile")==0 ){
00490
00491 if( dname != NULL )
00492 ERREX("*** Can't use 2 -dfile/-1Dfile options!") ;
00493 if( matvec )
00494 ERREX("*** Can't combine -dfile/-1Dfile with -matvec!") ;
00495 if( dcode > 0 || rotarg > 0 )
00496 ERREX("*** Can't use -dfile/-1Dfile with -shift or -rotate options!") ;
00497 if( rotpar_dset != NULL )
00498 ERREX("*** Can't use -dfile/-1Dfile with -rotparent!") ;
00499
00500 if( strcmp(argv[iopt],"-dfile") ==0 ) dmode = MODE_DFILE ;
00501 else if( strcmp(argv[iopt],"-1Dfile")==0 ) dmode = MODE_1DFILE ;
00502
00503 dname = argv[++iopt] ;
00504
00505 dim = get_dfile_params( dname , dmode ) ;
00506 if( dim == NULL )
00507 ERREX("*** Error with -dfile or -1Dfile!") ;
00508
00509 dar = MRI_FLOAT_PTR(dim) ;
00510 ndar = dim->ny ;
00511 iopt++ ; continue ;
00512 }
00513
00514 fprintf(stderr,"*** Unknown option: %s\n",argv[iopt]) ; exit(1) ;
00515 }
00516
00517
00518
00519 if( gridpar_dset != NULL && rotpar_dset == NULL ){
00520 fprintf(stderr,"+++ WARNING: -gridparent means nothing without -rotparent!\n");
00521 DSET_delete( gridpar_dset ) ;
00522 gridpar_dset = NULL ;
00523 }
00524
00525 if( gridpar_dset != NULL && dopoints ){
00526 fprintf(stderr,"+++ WARNING: -gridparent means nothing with -points!\n") ;
00527 DSET_delete( gridpar_dset ) ;
00528 gridpar_dset = NULL ;
00529 }
00530
00531 if( doorigin && !dopoints ){
00532 fprintf(stderr,"+++ WARNING: -origin means nothing without -points!\n") ;
00533 }
00534
00535 if( dar != NULL && dopoints )
00536 ERREX("*** Can't combine -dfile/-1Dfile with -points!") ;
00537
00538 if( matvec==0 && dcode<0 && rotarg<0 && rotpar_dset==NULL && dar==NULL )
00539 ERREX("Don't you want to do anything [no -rotate,-shift,-matvec,-rotparent,-dfile]?");
00540
00541
00542
00543 if( iopt >= argc && !dopoints ) ERREX("*** No input dataset?") ;
00544
00545 if( iopt < argc ){
00546 dset = THD_open_dataset( argv[iopt] ) ;
00547 if( dset == NULL ){
00548 fprintf(stderr,"*** Cannot open dataset %s!\n",argv[iopt]); exit(1);
00549 }
00550 } else {
00551 dset = EDIT_empty_copy(NULL) ;
00552 if( !doorigin )
00553 fprintf(stderr,"+++ WARNING: no -origin and no input dataset on command line!\n") ;
00554 }
00555
00556 if( dopoints && !doorigin ){
00557 xo = dset->daxes->xxorg + 0.5*(dset->daxes->nxx - 1)*dset->daxes->xxdel ;
00558 yo = dset->daxes->yyorg + 0.5*(dset->daxes->nyy - 1)*dset->daxes->yydel ;
00559 zo = dset->daxes->zzorg + 0.5*(dset->daxes->nzz - 1)*dset->daxes->zzdel ;
00560 fprintf(stderr,"+++ Using -origin %g %g %g\n",xo,yo,zo) ;
00561 }
00562
00563
00564
00565 ihand = THD_handedness(dset) ;
00566
00567 if( rotarg > 0 ){
00568 int neg ;
00569 iopt = rotarg ;
00570
00571 th1 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
00572 switch( *cpt ){
00573 default: fprintf(stderr,"*** Illegal code after th1 in -rotate\n");exit(1);
00574 case '\0': case 'x': case 'X': ax1 = 0 ; neg = 0 ; break ;
00575 case 'y': case 'Y': ax1 = 1 ; neg = 0 ; break ;
00576 case 'z': case 'Z': ax1 = 2 ; neg = 0 ; break ;
00577
00578 case 'A': case 'P':
00579 case 'R': case 'L':
00580 case 'I': case 'S': ax1 = THD_axcode(dset,*cpt) ;
00581 neg = (ax1 < 0) ;
00582 ax1 = abs(ax1) - 1 ; break ;
00583 }
00584 if( neg ) th1 = -th1 ;
00585
00586 th2 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
00587 switch( *cpt ){
00588 default: fprintf(stderr,"*** Illegal code after th2 in -rotate\n");exit(1);
00589 case 'x': case 'X': ax2 = 0 ; break ;
00590 case '\0': case 'y': case 'Y': ax2 = 1 ; break ;
00591 case 'z': case 'Z': ax2 = 2 ; break ;
00592
00593 case 'A': case 'P':
00594 case 'R': case 'L':
00595 case 'I': case 'S': ax2 = THD_axcode(dset,*cpt) ;
00596 neg = (ax2 < 0) ;
00597 ax2 = abs(ax2) - 1 ; break ;
00598 }
00599 if( neg ) th2 = -th2 ;
00600
00601 th3 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
00602 switch( *cpt ){
00603 default: fprintf(stderr,"*** Illegal code after th3 in -rotate\n");exit(1);
00604 case 'x': case 'X': ax3 = 0 ; break ;
00605 case 'y': case 'Y': ax3 = 1 ; break ;
00606 case '\0': case 'z': case 'Z': ax3 = 2 ; break ;
00607
00608 case 'A': case 'P':
00609 case 'R': case 'L':
00610 case 'I': case 'S': ax3 = THD_axcode(dset,*cpt) ;
00611 neg = (ax3 < 0) ;
00612 ax3 = abs(ax3) - 1 ; break ;
00613 }
00614 if( neg ) th3 = -th3 ;
00615
00616 if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ){
00617 fprintf(stderr,"*** Why are all the angles after -rotate equal to 0?\n") ;
00618 exit(1) ;
00619 }
00620
00621 if( ax1 < 0 || ax1 > 2 || ax2 < 0 || ax2 > 2 || ax3 < 0 || ax3 > 2 ){
00622 fprintf(stderr,"*** Some error occured: can't understand axes codes!\n") ;
00623 exit(1) ;
00624 }
00625
00626 if( ihand < 0 ){ th1 = -th1 ; th2 = -th2 ; th3 = -th3 ; }
00627
00628 #if 0
00629 fprintf(stderr,"ihand=%d th1=%g th2=%g th3=%g\n",ihand,th1,th2,th3);
00630 fprintf(stderr,"ax1=%d ax2=%d ax3=%d\n",ax1,ax2,ax3) ;
00631 #endif
00632
00633
00634
00635 if( dopoints )
00636 rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ;
00637 }
00638
00639
00640
00641 if( dcode > 0 && (cdx != '\0' || cdy != '\0' || cdz != '\0') ){
00642 float qdx=0 , qdy=0 , qdz=0 ;
00643
00644 adx = THD_axcode(dset,cdx) ;
00645 if( adx > -99 || dx != 0.0 ){
00646 switch( adx ){
00647 case 1: qdx = -dx ; break ;
00648 default:
00649 case -1: qdx = dx ; break ;
00650 case 2: qdy = -dx ; break ;
00651 case -2: qdy = dx ; break ;
00652 case 3: qdz = -dx ; break ;
00653 case -3: qdz = dx ; break ;
00654 }
00655 }
00656
00657 ady = THD_axcode(dset,cdy) ;
00658 if( ady > -99 || dy != 0.0 ){
00659 switch( ady ){
00660 case 1: qdx = -dy ; break ;
00661 case -1: qdx = dy ; break ;
00662 case 2: qdy = -dy ; break ;
00663 default:
00664 case -2: qdy = dy ; break ;
00665 case 3: qdz = -dy ; break ;
00666 case -3: qdz = dy ; break ;
00667 }
00668 }
00669
00670 adz = THD_axcode(dset,cdz) ;
00671 if( adz > -99 || dz != 0.0 ){
00672 switch( adz ){
00673 case 1: qdx = -dz ; break ;
00674 case -1: qdx = dz ; break ;
00675 case 2: qdy = -dz ; break ;
00676 case -2: qdy = dz ; break ;
00677 case 3: qdz = -dz ; break ;
00678 default:
00679 case -3: qdz = dz ; break ;
00680 }
00681 }
00682
00683 if( verb )
00684 fprintf(stderr,"Shifting parameters:\n"
00685 " direction codes: adx=%d ady=%d adz=%d\n"
00686 " input values: dx=%g dy=%g dz=%g\n"
00687 " output values: qdx=%g qdy=%g qdz=%g\n" ,
00688 adx,ady,adz , dx,dy,dz , qdx,qdy,qdz ) ;
00689
00690 dx = qdx ; dy = qdy ; dz = qdz ;
00691 }
00692
00693
00694
00695 if( matvec == MATVEC_DICOM ){
00696 pp = DBLE_mat_to_dicomm( dset ) ;
00697 ppt = TRANSPOSE_DMAT(pp);
00698 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp);
00699 tvec = DMATVEC(ppt,tvec);
00700 }
00701
00702
00703
00704
00705 if( gridpar_dset != NULL ){
00706 int mm , nz_gp , nz_ds ;
00707
00708
00709
00710 mm = THD_dataset_mismatch( gridpar_dset , dset ) ;
00711 if( mm & (MISMATCH_DELTA | MISMATCH_ORIENT) ){
00712 fprintf(stderr,"*** Fatal Error:\n"
00713 "*** -gridparent dataset and input dataset don't\n"
00714 "*** match in grid spacing and/or orientation!\n" ) ;
00715 exit(1) ;
00716 }
00717
00718 if( DSET_NX(gridpar_dset) != DSET_NX(dset) ||
00719 DSET_NY(gridpar_dset) != DSET_NY(dset) ){
00720
00721 fprintf(stderr,"*** Fatal Error:\n"
00722 "*** -gridparent and input datasets\n"
00723 "*** don't match in x,y dimensions!\n" ) ;
00724 exit(1) ;
00725 }
00726
00727
00728
00729 nz_gp = DSET_NZ(gridpar_dset) ; nz_ds = DSET_NZ(dset) ;
00730
00731 if( nz_gp < nz_ds ){
00732 fprintf(stderr,"*** Fatal Error:\n"
00733 "*** -gridparent has fewer slices than input dataset!\n") ;
00734 exit(1) ;
00735 }
00736
00737 if( nz_gp > nz_ds ){
00738 int npad1 = (nz_gp - nz_ds) / 2 ;
00739 int npad2 = (nz_gp - nz_ds) - npad1 ;
00740 int add_I=0, add_S=0, add_A=0, add_P=0, add_L=0, add_R=0 ;
00741 THD_3dim_dataset * pset ;
00742 char *sp1,*sp2 ;
00743
00744
00745
00746 switch( dset->daxes->zzorient ){
00747 case ORI_R2L_TYPE:
00748 case ORI_L2R_TYPE: add_R=npad1; add_L=npad2; sp1="R"; sp2="L"; break;
00749
00750 case ORI_P2A_TYPE:
00751 case ORI_A2P_TYPE: add_A=npad1; add_P=npad2; sp1="A"; sp2="P"; break;
00752
00753 case ORI_I2S_TYPE:
00754 case ORI_S2I_TYPE: add_I=npad1; add_S=npad2; sp1="I"; sp2="S"; break;
00755 }
00756
00757
00758
00759 if( verb )
00760 fprintf(stderr,"+++ Zero padding to match -gridparent: -%s %d -%s %d\n",
00761 sp1,npad1,sp2,npad2 ) ;
00762
00763 pset = THD_zeropad( dset,
00764 add_I,add_S,add_A,add_P,add_L,add_R,
00765 "Elvis" , ZPAD_PURGE ) ;
00766
00767 if( pset == NULL ){
00768 fprintf(stderr,"*** Fatal Error:\n"
00769 "*** Can't properly zeropad input dataset!\n" ) ;
00770 exit(1) ;
00771 }
00772
00773
00774
00775 DSET_delete(dset) ; dset = pset ;
00776 }
00777 }
00778
00779 if( rotpar_dset != NULL ){
00780 ATR_float * atr ;
00781 float * matar , sum ;
00782 THD_fvec3 fv ;
00783 THD_dfvec3 dv,ev,qv , cv_e2, cv_e1, cv_s1, cv_s2 ;
00784
00785
00786
00787 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ;
00788 matar = atr->fl ;
00789 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
00790 matar[4],matar[5],matar[6],
00791 matar[8],matar[9],matar[10] ) ;
00792 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
00793
00794
00795
00796 pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ;
00797 sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0])
00798 +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1])
00799 +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0);
00800 if( sum > 0.01 ) ERREX("-rotparent matrix not orthogonal!") ;
00801
00802
00803
00804
00805
00806
00807 fv = THD_dataset_center( dset ) ;
00808 FVEC3_TO_DFVEC3( fv , cv_e2 ) ;
00809
00810
00811
00812 if( gridpar_dset != NULL ){
00813 fv = THD_dataset_center( gridpar_dset ) ;
00814 FVEC3_TO_DFVEC3( fv , cv_e1 ) ;
00815 } else {
00816 cv_e1 = cv_e2 ;
00817 }
00818
00819
00820
00821 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_OLD" ) ;
00822 LOAD_DFVEC3( cv_s2 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
00823
00824
00825
00826 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_BASE" ) ;
00827 LOAD_DFVEC3( cv_s1 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
00828
00829
00830
00831
00832
00833 dv = SUB_DFVEC3( cv_e2 , cv_s2 ) ;
00834 ev = DMATVEC( rmat , dv ) ;
00835
00836 dv = ev ;
00837
00838 ev = SUB_DFVEC3( cv_e1 , cv_s1 ) ;
00839
00840 qv = SUB_DFVEC3( dv , ev ) ;
00841
00842 tvec = ADD_DFVEC3( tvec , qv ) ;
00843
00844
00845
00846 pp = DBLE_mat_to_dicomm( dset ) ;
00847 ppt = TRANSPOSE_DMAT(pp);
00848 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp); tvec = DMATVEC(ppt,tvec);
00849
00850
00851
00852 if( gridpar_dset != NULL ){
00853 dset->daxes->xxorg = gridpar_dset->daxes->xxorg ;
00854 dset->daxes->yyorg = gridpar_dset->daxes->yyorg ;
00855 dset->daxes->zzorg = gridpar_dset->daxes->zzorg ;
00856
00857
00858
00859 if( dset->taxis != NULL && dset->taxis->nsl > 0 ){
00860 dset->taxis->zorg_sl = dset->daxes->zzorg ;
00861 }
00862 }
00863
00864 matvec = MATVEC_ORDER ;
00865 }
00866
00867
00868
00869 if( dopoints ){
00870 THD_dfvec3 xyzorg ;
00871 if( !matvec ) LOAD_DFVEC3( tvec , dx,dy,dz ) ;
00872 if( dcode < 0 ) dcode = DELTA_AFTER ;
00873 LOAD_DFVEC3(xyzorg,xo,yo,zo) ;
00874 rotate_stdin_points( xyzorg , rmat , dcode,tvec ) ;
00875 exit(0) ;
00876 }
00877
00878
00879
00880 if( dset->taxis != NULL && dset->taxis->nsl > 0 ){
00881 int ndz ;
00882 int kk,jj , nsl = dset->taxis->nsl ;
00883
00884 if( matvec )
00885 ndz = (int) rint( tvec.xyz[2] / fabs(dset->daxes->zzdel) ) ;
00886 else
00887 ndz = (int) rint( dz / fabs(dset->daxes->xxdel) ) ;
00888
00889 if( ndz != 0 ){
00890 float * tsl = (float *)malloc(sizeof(float)*nsl) ;
00891 for( kk=0 ; kk < nsl ; kk ++ ){
00892 jj = kk - ndz ;
00893 if( jj < 0 || jj >= nsl ) tsl[kk] = 0.0 ;
00894 else tsl[kk] = dset->taxis->toff_sl[jj] ;
00895 }
00896 EDIT_dset_items( dset , ADN_toff_sl , tsl , ADN_none ) ;
00897 free(tsl) ;
00898 if( verb )
00899 fprintf(stderr,"+++ adjusting time-offsets by %d slices\n",ndz) ;
00900 }
00901 }
00902
00903
00904
00905 DSET_mallocize(dset) ;
00906 DSET_load(dset) ;
00907 dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ;
00908
00909 dset->idcode = MCW_new_idcode() ;
00910 EDIT_dset_items( dset ,
00911 ADN_prefix , new_prefix ,
00912 ADN_label1 , new_prefix ,
00913 ADN_none ) ;
00914 if( THD_is_file(dset->dblk->diskptr->header_name) ){
00915 fprintf(stderr,
00916 "** ERROR: Output file %s already exists -- cannot continue!\n",
00917 dset->dblk->diskptr->header_name ) ;
00918 exit(1) ;
00919 }
00920
00921
00922
00923 tross_Make_History( "3drotate" , argc,argv , dset ) ;
00924
00925 nvox = DSET_NVOX(dset) ;
00926 fvol = (float *) malloc( sizeof(float) * nvox ) ;
00927
00928 nval = DSET_NVALS(dset) ;
00929 if( verb ){
00930 fprintf(stderr,"+++ %d sub-bricks: ",nval) ;
00931 cputim = COX_cpu_time() ;
00932 }
00933
00934
00935
00936 { THD_fvec3 cv ;
00937 cv = THD_dataset_center( dset ) ;
00938 THD_set_float_atr( dset->dblk , "ROTATE_CENTER_OLD" , 3 , cv.xyz ) ;
00939 THD_set_float_atr( dset->dblk , "ROTATE_CENTER_BASE" , 3 , cv.xyz ) ;
00940 }
00941
00942
00943
00944 for( ival=0 ; ival < nval ; ival++ ){
00945
00946 if( verb ) fprintf(stderr,"%d",ival) ;
00947
00948 EDIT_coerce_type( nvox ,
00949 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ,
00950 MRI_float,fvol ) ;
00951
00952 if( verb ) fprintf(stderr,".") ;
00953
00954 if( clipit ){
00955 register int ii ; register float bb,tt ;
00956 bb = tt = fvol[0] ;
00957 for( ii=1 ; ii < nvox ; ii++ ){
00958 if( fvol[ii] < bb ) bb = fvol[ii] ;
00959 else if( fvol[ii] > tt ) tt = fvol[ii] ;
00960 }
00961 cbot = bb ; ctop = tt ;
00962 }
00963
00964
00965
00966 skipit = 0 ;
00967 if( dar != NULL ){
00968 THD_dvecmat dvm ; char rotcom[256] ; int jj=ival ;
00969 static int ndar_over=0 ;
00970
00971 if( jj >= ndar ){
00972 jj = ndar-1 ;
00973 if( ndar_over == 0 )
00974 fprintf(stderr,
00975 "++ WARNING: from brick %d on, using last line (%d) from %s\n",
00976 jj , ndar-1 , dname ) ;
00977 ndar_over++ ;
00978 }
00979
00980 sprintf(rotcom,"-rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP",
00981 dar[0+6*jj] , dar[1+6*jj] , dar[2+6*jj] ,
00982 dar[3+6*jj] , dar[4+6*jj] , dar[5+6*jj] ) ;
00983
00984 dvm = THD_rotcom_to_matvec( dset , rotcom ) ;
00985
00986 rmat = dvm.mm ; tvec = dvm.vv ; matvec = MATVEC_ORDER ;
00987
00988 skipit = (dar[0+6*jj]==0.0 && dar[1+6*jj]==0.0 && dar[2+6*jj]==0.0 &&
00989 dar[3+6*jj]==0.0 && dar[4+6*jj]==0.0 && dar[5+6*jj]==0.0 );
00990
00991 if( !skipit ){
00992 skipit = ( fabs(rmat.mat[0][0]-1.0) < 0.00001 ) &&
00993 ( fabs(rmat.mat[1][1]-1.0) < 0.00001 ) &&
00994 ( fabs(rmat.mat[2][2]-1.0) < 0.00001 ) &&
00995 ( fabs(tvec.xyz[0]) < 0.001 ) &&
00996 ( fabs(tvec.xyz[1]) < 0.001 ) &&
00997 ( fabs(tvec.xyz[2]) < 0.001 ) ;
00998 }
00999
01000 if( verb && skipit ) fprintf(stderr,"[skip]");
01001 }
01002
01003
01004
01005 if( !skipit ){
01006 if( matvec ){
01007 THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
01008 fabs(DSET_DX(dset)) ,
01009 fabs(DSET_DY(dset)) ,
01010 fabs(DSET_DZ(dset)) ,
01011 fvol , rmat , tvec ) ;
01012 } else {
01013 THD_rota_vol( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
01014 fabs(DSET_DX(dset)),
01015 fabs(DSET_DY(dset)),
01016 fabs(DSET_DZ(dset)), fvol ,
01017 ax1,th1, ax2,th2, ax3,th3, dcode,dx,dy,dz ) ;
01018
01019
01020
01021 rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ;
01022 LOAD_DFVEC3( tvec , dx,dy,dz ) ;
01023 if( dcode == DELTA_BEFORE ) tvec = DMATVEC(rmat,tvec) ;
01024 }
01025 }
01026
01027
01028
01029 { float matar[12] ;
01030 THD_dmat33 drmat ;
01031 THD_dfvec3 dvec ;
01032 char anam[64] ;
01033
01034
01035
01036
01037 pp = DBLE_mat_to_dicomm(dset) ; ppt = TRANSPOSE_DMAT(pp) ;
01038 drmat = DMAT_MUL(pp,rmat) ; drmat = DMAT_MUL(drmat,ppt);
01039 dvec = DMATVEC (pp,tvec) ;
01040
01041 UNLOAD_DMAT(drmat,matar[0],matar[1],matar[2],
01042 matar[4],matar[5],matar[6],
01043 matar[8],matar[9],matar[10] ) ;
01044 UNLOAD_DFVEC3(dvec,matar[3],matar[7],matar[11]) ;
01045 sprintf(anam,"ROTATE_MATVEC_%06d",ival) ;
01046 THD_set_float_atr( dset->dblk , anam , 12 , matar ) ;
01047 }
01048
01049 if( verb ) fprintf(stderr,".") ;
01050
01051 if( clipit ){
01052 register int ii ; register float bb,tt ;
01053 bb = cbot ; tt = ctop ;
01054 for( ii=0 ; ii < nvox ; ii++ ){
01055 if( fvol[ii] < bb ) fvol[ii] = bb ;
01056 else if( fvol[ii] > tt ) fvol[ii] = tt ;
01057 }
01058 }
01059
01060 EDIT_coerce_type( nvox , MRI_float,fvol ,
01061 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ) ;
01062
01063 }
01064
01065 if( verb ){
01066 cputim = COX_cpu_time() - cputim ;
01067 fprintf(stderr,"\n+++ CPU time=%10.3g s" , cputim) ;
01068 if( nval > 1 ) fprintf(stderr," [= %10.3g s/sub-brick]" , cputim/nval) ;
01069 fprintf(stderr,"\n+++ Writing dataset to disk in %s", DSET_BRIKNAME(dset) ) ;
01070 }
01071
01072 dset->dblk->master_nvals = 0 ;
01073 DSET_write(dset) ;
01074 if( verb ) fprintf(stderr,"\n") ;
01075 exit(0) ;
01076 }
01077
01078
01079
01080 #include <ctype.h>
01081
01082 #define NBUF 1024
01083
01084 #define DUPOUT(n) fprintf(fpout,"%s",linbuf+n)
01085
01086 static void rotate_stdin_points( THD_dfvec3 xyzorg, THD_dmat33 rmat,
01087 int dcode, THD_dfvec3 tvec )
01088 {
01089 char linbuf[NBUF] , *cp ;
01090 FILE *fpin=stdin , *fpout=stdout ;
01091 int ii , kk , nbuf , ll , nn , ld ;
01092 double xx,yy,zz ;
01093 THD_dfvec3 xyz ;
01094
01095 if( verb ){
01096 DUMP_DMAT33("Rotation",rmat) ;
01097 }
01098
01099
01100
01101
01102
01103 ll = ld = 0 ;
01104 while(1){
01105 ll++ ;
01106 cp = fgets( linbuf , NBUF , fpin ) ;
01107 if( cp == NULL ) break ;
01108 kk = strlen(linbuf) ;
01109 if( kk == 0 ) continue ;
01110
01111
01112
01113 for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ;
01114 if( ii == kk ||
01115 (linbuf[ii] == '/' && linbuf[ii+1] == '/') ){
01116
01117 DUPOUT(0) ; continue ;
01118 }
01119
01120
01121
01122 nn = sscanf(linbuf+ii , "%lf%lf%lf%n" , &xx,&yy,&zz,&nbuf ) ;
01123 if( nn < 3 ){
01124 fprintf(stderr,"+++ WARNING: input line %d was incomplete\n",ll) ;
01125 continue ;
01126 }
01127 nbuf += ii ;
01128
01129
01130
01131 LOAD_DFVEC3(xyz , xx,yy,zz) ;
01132 xyz = SUB_DFVEC3(xyz,xyzorg) ;
01133 if( dcode == DELTA_BEFORE ) xyz = ADD_DFVEC3(xyz,tvec) ;
01134 xyz = DMATVEC(rmat,xyz) ;
01135 if( dcode == DELTA_AFTER ) xyz = ADD_DFVEC3(xyz,tvec) ;
01136 xyz = ADD_DFVEC3(xyz,xyzorg) ;
01137
01138 fprintf(fpout,"%g %g %g%s",xyz.xyz[0],xyz.xyz[1],xyz.xyz[2],linbuf+nbuf) ;
01139 ld++ ;
01140
01141 }
01142
01143 if( verb )
01144 fprintf(stderr,"-points: read %d lines, wrote %d lines\n",ll-1,ld) ;
01145 }
01146
01147
01148
01149
01150
01151
01152 MRI_IMAGE * get_dfile_params( char * fname , int mode )
01153 {
01154 MRI_IMAGE *outim , *flim ;
01155 float *oar , *far ;
01156 int nx,ny , ii,jj ;
01157
01158 if( fname == NULL || fname[0] == '\0' ) return NULL ;
01159
01160 flim = mri_read_ascii( fname ) ;
01161 if( flim == NULL ) return NULL ;
01162
01163 nx = flim->nx ; ny = flim->ny ; far = MRI_FLOAT_PTR(flim) ;
01164
01165 if( mode == MODE_DFILE ){
01166 if( nx < 7 ){
01167 fprintf(stderr,"** -dfile %s has too few columns!\n",fname) ;
01168 mri_free(flim) ; return NULL ;
01169 }
01170 outim = mri_new( 6 , ny , MRI_float ) ;
01171 oar = MRI_FLOAT_PTR(outim) ;
01172
01173 for( jj=0 ; jj < ny ; jj++ ){
01174 oar[0+6*jj] = far[1+nx*jj] ;
01175 oar[1+6*jj] = far[2+nx*jj] ;
01176 oar[2+6*jj] = far[3+nx*jj] ;
01177 oar[3+6*jj] = far[4+nx*jj] ;
01178 oar[4+6*jj] = far[5+nx*jj] ;
01179 oar[5+6*jj] = far[6+nx*jj] ;
01180 }
01181
01182 mri_free(flim) ;
01183
01184 } else if( mode == MODE_1DFILE ){
01185
01186 if( nx < 6 ){
01187 fprintf(stderr,"** -1Dfile %s has too few columns!\n",fname) ;
01188 mri_free(flim) ; return NULL ;
01189 } else if( nx == 6 ){
01190 outim = flim ;
01191 } else {
01192 outim = mri_new( 6 , ny , MRI_float ) ;
01193 oar = MRI_FLOAT_PTR(outim) ;
01194 for( jj=0 ; jj < ny ; jj++ ){
01195 oar[0+6*jj] = far[0+nx*jj] ;
01196 oar[1+6*jj] = far[1+nx*jj] ;
01197 oar[2+6*jj] = far[2+nx*jj] ;
01198 oar[3+6*jj] = far[3+nx*jj] ;
01199 oar[4+6*jj] = far[4+nx*jj] ;
01200 oar[5+6*jj] = far[5+nx*jj] ;
01201 }
01202 mri_free(flim) ;
01203 }
01204
01205 } else {
01206 fprintf(stderr,"** get_dfile_params: illegal mode=%d\n",mode) ;
01207 mri_free(flim) ; return NULL ;
01208 }
01209
01210 return outim ;
01211 }
