Doxygen Source Code Documentation
3drotate.c File Reference
#include "thd_shear3d.h"
#include <ctype.h>
Go to the source code of this file.
Defines | |
#define | MATVEC_DICOM 1 |
#define | MATVEC_ORDER 2 |
#define | MODE_DFILE 1 |
#define | MODE_1DFILE 2 |
#define | ERREX(str) (fprintf(stderr,"*** %s\n",str),exit(1)) |
#define | NBUF 1024 |
#define | DUPOUT(n) fprintf(fpout,"%s",linbuf+n) |
Functions | |
void | rotate_stdin_points (THD_dfvec3, THD_dmat33, int, THD_dfvec3) |
MRI_IMAGE * | get_dfile_params (char *fname, int mode) |
int | main (int argc, char *argv[]) |
Variables | |
int | verb = 0 |
Define Documentation
|
Definition at line 1084 of file 3drotate.c. Referenced by rotate_stdin_points(). |
|
|
|
Definition at line 15 of file 3drotate.c. Referenced by main(). |
|
Definition at line 16 of file 3drotate.c. Referenced by main(). |
|
Definition at line 23 of file 3drotate.c. Referenced by get_dfile_params(), and main(). |
|
Definition at line 22 of file 3drotate.c. Referenced by get_dfile_params(), and main(). |
|
02 May 2005: record operation in header * Definition at line 1082 of file 3drotate.c. Referenced by rotate_stdin_points(). |
Function Documentation
|
rmat = DMAT_INV(rmat) ; * Definition at line 1152 of file 3drotate.c. References far, MODE_1DFILE, MODE_DFILE, MRI_FLOAT_PTR, mri_free(), mri_new(), mri_read_ascii(), MRI_IMAGE::nx, and MRI_IMAGE::ny. Referenced by main().
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 ){ /* skip 1st element of each row */ 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 ){ /* first 6 elements of each row */ 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 } |
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 28 of file 3drotate.c. References abs, ADD_DFVEC3, addto_args(), AFNI_logger(), argc, THD_3dim_dataset::daxes, DBLE_mat_to_dicomm(), THD_3dim_dataset::dblk, DELTA_AFTER, DELTA_BEFORE, DELTA_FIXED, DMAT_MUL, DMATVEC, DSET_delete, DSET_NX, DSET_NY, DSET_NZ, EDIT_empty_copy(), ev, ATR_float::fl, FVEC3_TO_DFVEC3, get_dfile_params(), LOAD_DFVEC3, LOAD_DIAG_DMAT, LOAD_DMAT, machdep(), mainENTRY, MASTER_SHORTHELP_STRING, THD_dmat33::mat, MATVEC_DICOM, MATVEC_ORDER, MISMATCH_DELTA, MISMATCH_ORIENT, MODE_1DFILE, MODE_DFILE, MRI_CUBIC, MRI_FLOAT_PTR, MRI_FOURIER, MRI_FOURIER_NOPAD, mri_free(), MRI_HEPTIC, MRI_LINEAR, MRI_NN, MRI_QUINTIC, mri_read_ascii(), ATR_float::nfl, THD_timeaxis::nsl, MRI_IMAGE::nx, THD_dataxes::nxx, MRI_IMAGE::ny, THD_dataxes::nyy, THD_dataxes::nzz, ORI_A2P_TYPE, ORI_I2S_TYPE, ORI_L2R_TYPE, ORI_P2A_TYPE, ORI_R2L_TYPE, ORI_S2I_TYPE, PRINT_VERSION, rot_to_matrix(), strtod(), SUB_DFVEC3, THD_3dim_dataset::taxis, THD_axcode(), THD_dataset_center(), THD_dataset_mismatch(), THD_find_float_atr(), THD_handedness(), THD_open_dataset(), THD_open_one_dataset(), THD_rota_method(), THD_rota_setpad(), THD_zeropad(), TRANSPOSE_DMAT, verb, THD_dataxes::xxdel, THD_dataxes::xxorg, THD_dataxes::yydel, THD_dataxes::yyorg, THD_timeaxis::zorg_sl, ZPAD_PURGE, THD_dataxes::zzdel, THD_dataxes::zzorg, and THD_dataxes::zzorient.
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 ; /* 11 Apr 2000 and 16 Apr 2002 */ 00040 float cbot,ctop ; 00041 00042 int matvec=0 ; /* 19 July 2000 */ 00043 THD_dmat33 rmat , pp,ppt ; 00044 THD_dfvec3 tvec ; 00045 00046 int dopoints=0 , doorigin=0 ; /* 21 Nov 2000 */ 00047 double xo=0.0 , yo=0.0 , zo=0.0 ; 00048 00049 THD_3dim_dataset *rotpar_dset=NULL , *gridpar_dset=NULL ; /* 07 Feb 2001 */ 00050 00051 int zpad=0 ; /* 05 Feb 2001 */ 00052 00053 char *dname=NULL ; /* 19 Jun 2001 */ 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 /*-- read command line arguments --*/ 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 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/ 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 ){ /* 05 Feb 2001 */ 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 ){ /* 21 Nov 2000 */ 00278 dopoints = 1 ; 00279 iopt++ ; continue ; 00280 } 00281 00282 if( strncmp(argv[iopt],"-origin",7) == 0 ){ /* 21 Nov 2000 */ 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 ){ /* 07 Feb 2001 */ 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 ){ /* 07 Feb 2001 */ 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 ){ /* 19 Jul 2000 */ 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){ /* 20 July 2000 */ 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 /* check if matrix is approximately orthogonal */ 00371 /* [will be orthogonalized in rot_to_shear_matvec() in thd_shear3d.c] */ 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 ){ /* 11 Apr 2000 */ 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 ){ /* 16 Apr 2002 */ 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 ){ /* 13 May 2003 */ 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 ; /* save and process later */ 00484 iopt += 4 ; continue ; 00485 } 00486 00487 /*-- 19 Jun 2001 --*/ 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) ; /* 6 x ndar array */ 00510 ndar = dim->ny ; 00511 iopt++ ; continue ; 00512 } 00513 00514 fprintf(stderr,"*** Unknown option: %s\n",argv[iopt]) ; exit(1) ; 00515 } 00516 00517 /*-- check for legal combinations --*/ 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 /** read input dataset */ 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) ; /* 21 Nov 2000: need a fake dataset */ 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 /* now can process rotation arguments */ 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 /* notice the minus signs on the angles ------+ */ 00634 /* | | | */ 00635 if( dopoints ) /* V V V */ 00636 rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ; /* 21 Nov 2000 */ 00637 } 00638 00639 /* may need to process shift arguments as well */ 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 ){ /* 29 Jan 2001: skip if purely 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 ){ /* 29 Jan 2001 */ 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 ){ /* 29 Jan 2001 */ 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 /*- 19 July 2000: now can deal with -matvec_dicom case -*/ 00694 00695 if( matvec == MATVEC_DICOM ){ /* convert matrix/vector */ 00696 pp = DBLE_mat_to_dicomm( dset ) ; /* to dataset coord order */ 00697 ppt = TRANSPOSE_DMAT(pp); /* from the DICOM order! */ 00698 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp); 00699 tvec = DMATVEC(ppt,tvec); 00700 } 00701 00702 /*-- 07 Feb 2001: deal with -rotparent and -gridparent case 00703 [very similar to -matvec_dicom, actually] --*/ 00704 00705 if( gridpar_dset != NULL ){ 00706 int mm , nz_gp , nz_ds ; 00707 00708 /* check for compatibility! */ 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 /* check for zero padding */ 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 /* where to add slices? and how many? */ 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 /* add them on */ 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 /* toss input datset, replace with padded one */ 00774 00775 DSET_delete(dset) ; dset = pset ; 00776 } 00777 } 00778 00779 if( rotpar_dset != NULL ){ /* compute transformation from -rotparent */ 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 /* load transformation from rotparent */ 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 /* check if matrix is orthogonal */ 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 /* must alter shift [tvec] to allow for differing 00803 coordinates in the rotparent, gridparent, and input datasets */ 00804 00805 /* cv_e2 = center of input dataset [Dicom coordinates] */ 00806 00807 fv = THD_dataset_center( dset ) ; /* dataset coords */ 00808 FVEC3_TO_DFVEC3( fv , cv_e2 ) ; /* convert to double */ 00809 00810 /* cv_e1 = center of gridparent */ 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 ; /* what else to do? */ 00817 } 00818 00819 /* cv_s2 = center of rotation in rotparent */ 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 /* cv_s1 = center of base dataset for rotparent */ 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 /* compute extra shift due to difference in 00830 center of rotation between rotparent and input dataset, 00831 then add in shifts caused by -twodup for rotparent and input */ 00832 00833 dv = SUB_DFVEC3( cv_e2 , cv_s2 ) ; 00834 ev = DMATVEC( rmat , dv ) ; /* R[E2-S2] */ 00835 00836 dv = ev ; /* vestige of a stupid bug, since fixed */ 00837 00838 ev = SUB_DFVEC3( cv_e1 , cv_s1 ) ; /* E1-S1 */ 00839 00840 qv = SUB_DFVEC3( dv , ev ) ; /* R[E2-S2] + S1-E1 */ 00841 00842 tvec = ADD_DFVEC3( tvec , qv ) ; /* shifted translation vector */ 00843 00844 /* convert transformation from Dicom to dataset coords */ 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 /* modify origin of output dataset to match -gridparent */ 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 /* 12 Feb 2001: adjust origin of time-offsets as well */ 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 ; /* flag that transform comes from rmat/tvec */ 00865 } 00866 00867 /*-- 21 Nov 2000: read (x,y,z) points from stdin, process them, quit --*/ 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 ) ; /* at end of file */ 00875 exit(0) ; 00876 } 00877 00878 /*-- 12 Feb 2001: adjust time-offsets for slice direction shifts --*/ 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) ) ; /* shift */ 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 /*- read dataset, prepare to process it, write back out (with new name) -*/ 00904 00905 DSET_mallocize(dset) ; 00906 DSET_load(dset) ; 00907 dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */ 00908 00909 dset->idcode = MCW_new_idcode() ; /* 08 Jun 1999 - is a new dataset */ 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 /* old history is already in the dataset */ 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 /* 03 May 2005: save rotation center */ 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 /*-- loop over all sub-bricks: copy into fvol, rotate fvol, copy back --*/ 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 ){ /* 11 Apr 2000 */ 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 /* 19 Jun 2001: get matrix/vector from rot params in file */ 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 ) ; /* thd_rotangles.c */ 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 /** carry out the rotation **/ 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 /* 02 May 2005: save matrix+vector for recording below */ 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 /** 02 May 2005: record operation in header **/ 01028 01029 { float matar[12] ; 01030 THD_dmat33 drmat ; 01031 THD_dfvec3 dvec ; 01032 char anam[64] ; 01033 01034 /* must transform [rmat,tvec] from dataset to DICOM coords */ 01035 /* (the inverse of the DICOM to dataset transforms earlier) */ 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 ){ /* 11 Apr 2000 */ 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 } /* end of loop over sub-bricks */ 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 ; /* 11 Apr 2000 hack */ 01073 DSET_write(dset) ; 01074 if( verb ) fprintf(stderr,"\n") ; 01075 exit(0) ; 01076 } |
|
Definition at line 1086 of file 3drotate.c. References ADD_DFVEC3, DELTA_AFTER, DELTA_BEFORE, DMATVEC, DUMP_DMAT33, DUPOUT, LOAD_DFVEC3, NBUF, SUB_DFVEC3, and THD_dfvec3::xyz.
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 /** rmat = DMAT_INV(rmat) ; **/ 01100 01101 /*-- loop over input lines --*/ 01102 01103 ll = ld = 0 ; 01104 while(1){ 01105 ll++ ; /* line count */ 01106 cp = fgets( linbuf , NBUF , fpin ) ; /* read the line */ 01107 if( cp == NULL ) break ; /* end of file => end of loop */ 01108 kk = strlen(linbuf) ; 01109 if( kk == 0 ) continue ; /* empty line => get next line */ 01110 01111 /* find 1st nonblank */ 01112 01113 for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ; /* nada */ 01114 if( ii == kk || /* all blanks */ 01115 (linbuf[ii] == '/' && linbuf[ii+1] == '/') ){ /* or comment */ 01116 01117 DUPOUT(0) ; continue ; 01118 } 01119 01120 /* scan line for data */ 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 ; /* position of next character after zz */ 01128 01129 /* process vector */ 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 } /* end of loop over input lines */ 01142 01143 if( verb ) 01144 fprintf(stderr,"-points: read %d lines, wrote %d lines\n",ll-1,ld) ; 01145 } |
Variable Documentation
|
Definition at line 18 of file 3drotate.c. Referenced by main(). |