/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2000, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ #include "thd_shear3d.h" /*-- prototypes for funcs at end of file --*/ static void rotate_stdin_points( THD_dfvec3, THD_dmat33, int,THD_dfvec3 ) ; /*------------------------------------------------------------------------------*/ #define MATVEC_DICOM 1 #define MATVEC_ORDER 2 static int verb=0 ; /*-- 19 Jun 2001 stuff --*/ #define MODE_DFILE 1 #define MODE_1DFILE 2 MRI_IMAGE * get_dfile_params( char * fname , int mode ) ; /*------------------------------------------------------------------------------*/ int main( int argc , char * argv[] ) { THD_3dim_dataset * dset=NULL ; char * new_prefix = "rota" , * cpt ; float dx=0 , dy=0 , dz=0 ; int ax1=0,ax2=1,ax3=2 , adx,ady,adz ; char cdx='\0',cdy='\0',cdz='\0' ; float th1=0.0,th2=0.0,th3=0.0 ; int iopt , nvox , rotarg=-1 , dcode=-1 , ival,nval , ihand ; float * fvol ; double cputim=0.0 ; int clipit=1 ; /* 11 Apr 2000 and 16 Apr 2002 */ float cbot=0.0,ctop=0.0 ; int matvec=0 ; /* 19 July 2000 */ THD_dmat33 rmat , pp,ppt ; THD_dfvec3 tvec ; int dopoints=0 , doorigin=0 ; /* 21 Nov 2000 */ double xo=0.0 , yo=0.0 , zo=0.0 ; THD_3dim_dataset *rotpar_dset=NULL , *gridpar_dset=NULL ; /* 07 Feb 2001 */ int zpad=0 ; /* 05 Feb 2001 */ char *dname=NULL ; /* 19 Jun 2001 */ int dmode=0 ; MRI_IMAGE *dim=NULL ; float *dar=NULL ; int ndar=0 ; int skipit=0 ; /*-------------------------------*/ LOAD_DIAG_DMAT(rmat,1.0,1.0,1.0) ; memset(&tvec, 0, sizeof(THD_dfvec3)); /*-- read command line arguments --*/ if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){ printf( "Usage: 3drotate [options] dataset\n" "Rotates and/or translates all bricks from an AFNI dataset.\n" "'dataset' may contain a sub-brick selector list.\n" "\n" "GENERIC OPTIONS:\n" " -prefix fname = Sets the output dataset prefix name to be 'fname'\n" " -verbose = Prints out progress reports (to stderr)\n" "\n" "OPTIONS TO SPECIFY THE ROTATION/TRANSLATION:\n" "-------------------------------------------\n" "*** METHOD 1 = direct specification:\n" "At most one of these shift options can be used:\n" " -ashift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n" " AFTER rotation.\n" " -bshift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n" " BEFORE rotation.\n" " The shift distances by default are along the (x,y,z) axes of the dataset\n" " storage directions (see the output of '3dinfo dataset'). To specify them\n" " anatomically, you can suffix a distance with one of the symbols\n" " 'R', 'L', 'A', 'P', 'I', and 'S', meaning 'Right', 'Left', 'Anterior',\n" " 'Posterior', 'Inferior', and 'Superior', respectively.\n" "\n" " -rotate th1 th2 th3\n" " Specifies the 3D rotation to be composed of 3 planar rotations:\n" " 1) 'th1' degrees about the 1st axis, followed by\n" " 2) 'th2' degrees about the (rotated) 2nd axis, followed by\n" " 3) 'th3' degrees about the (doubly rotated) 3rd axis.\n" " Which axes are used for these rotations is specified by placing\n" " one of the symbols 'R', 'L', 'A', 'P', 'I', and 'S' at the end\n" " of each angle (e.g., '10.7A'). These symbols denote rotation\n" " about the 'Right-to-Left', 'Left-to-Right', 'Anterior-to-Posterior',\n" " 'Posterior-to-Anterior', 'Inferior-to-Superior', and\n" " 'Superior-to-Inferior' axes, respectively. A positive rotation is\n" " defined by the right-hand rule.\n" "\n" "*** METHOD 2 = copy from output of 3dvolreg:\n" " -rotparent rset\n" " Specifies that the rotation and translation should be taken from the\n" " first 3dvolreg transformation found in the header of dataset 'rset'.\n" " -gridparent gset\n" " Specifies that the output dataset of 3drotate should be shifted to\n" " match the grid of dataset 'gset'. Can only be used with -rotparent.\n" " This dataset should be one this is properly aligned with 'rset' when\n" " overlaid in AFNI.\n" " * If -rotparent is used, then don't use -matvec, -rotate, or -[ab]shift.\n" " * If 'gset' has a different number of slices than the input dataset,\n" " then the output dataset will be zero-padded in the slice direction\n" " to match 'gset'.\n" " * These options are intended to be used to align datasets between sessions:\n" " S1 = SPGR from session 1 E1 = EPI from session 1\n" " S2 = SPGR from session 2 E2 = EPI from session 2\n" " 3dvolreg -twopass -twodup -base S1+orig -prefix S2reg S2+orig\n" " 3drotate -rotparent S2reg+orig -gridparent E1+orig -prefix E2reg E2+orig\n" " The result will have E2reg rotated from E2 in the same way that S2reg\n" " was from S2, and also shifted/padded (as needed) to overlap with E1.\n" "\n" "*** METHOD 3 = give the transformation matrix/vector directly:\n" " -matvec_dicom mfile\n" " -matvec_order mfile\n" " Specifies that the rotation and translation should be read from file\n" " 'mfile', which should be in the format\n" " u11 u12 u13 v1\n" " u21 u22 u23 v2\n" " u31 u32 u33 u3\n" " where each 'uij' and 'vi' is a number. The 3x3 matrix [uij] is the\n" " orthogonal matrix of the rotation, and the 3-vector [vi] is the -ashift\n" " vector of the translation.\n" "\n" "*** METHOD 4 = copy the transformation from 3dTagalign:\n" " -matvec_dset mset\n" " Specifies that the rotation and translation should be read from\n" " the .HEAD file of dataset 'mset', which was created by program\n" " 3dTagalign.\n" " * If -matvec_dicom is used, the matrix and vector are given in Dicom\n" " coordinate order (+x=L, +y=P, +z=S). This is the option to use\n" " if mfile is generated using 3dTagalign -matvec mfile.\n" " * If -matvec_order is used, the matrix and vector are given in the\n" " coordinate order of the dataset axes, whatever they may be.\n" " * You can't mix -matvec_* options with -rotate and -*shift.\n" "\n" "*** METHOD 5 = input rotation+shift parameters from an ASCII file:\n" " -dfile dname *OR* -1Dfile dname\n" " With these methods, the movement parameters for each sub-brick\n" " of the input dataset are read from the file 'dname'. This file\n" " should consist of columns of numbers in ASCII format. Six (6)\n" " numbers are read from each line of the input file. If the\n" " '-dfile' option is used, each line of the input should be at\n" " least 7 numbers, and be of the form\n" " ignored roll pitch yaw dS dL dP\n" " If the '-1Dfile' option is used, then each line of the input\n" " should be at least 6 numbers, and be of the form\n" " roll pitch yaw dS dL dP\n" " (These are the forms output by the '-dfile' and\n" " '-1Dfile' options of program 3dvolreg; see that\n" " program's -help output for the hideous details.)\n" " The n-th sub-brick of the input dataset will be transformed\n" " using the parameters from the n-th line of the dname file.\n" " If the dname file doesn't contain as many lines as the\n" " input dataset has sub-bricks, then the last dname line will\n" " be used for all subsequent sub-bricks. Excess columns or\n" " rows will be ignored.\n" " N.B.: Rotation is always about the center of the volume.\n" " If the parameters are derived from a 3dvolreg run\n" " on a dataset with a different center in xyz-space,\n" " the results may not be what you want!\n" " N.B.: You can't use -dfile/-1Dfile with -points (infra).\n" "\n" "POINTS OPTIONS (instead of datasets):\n" "------------------------------------\n" " -points\n" " -origin xo yo zo\n" " These options specify that instead of rotating a dataset, you will\n" " be rotating a set of (x,y,z) points. The points are read from stdin.\n" " * If -origin is given, the point (xo,yo,zo) is used as the center for\n" " the rotation.\n" " * If -origin is NOT given, and a dataset is given at the end of the\n" " command line, then the center of the dataset brick is used as\n" " (xo,yo,zo). The dataset will NOT be rotated if -points is given.\n" " * If -origin is NOT given, and NO dataset is given at the end of the\n" " command line, then xo=yo=zo=0 is assumed. You probably don't\n" " want this.\n" " * (x,y,z) points are read from stdin as 3 ASCII-formatted numbers per\n" " line, as in 3dUndump. Any succeeding numbers on input lines will\n" " be copied to the output, which will be written to stdout.\n" " * The input (x,y,z) coordinates are taken in the same order as the\n" " axes of the input dataset. If there is no input dataset, then\n" " negative x = R positive x = L }\n" " negative y = A positive y = P } e.g., the DICOM order\n" " negative z = I positive z = S }\n" " One way to dump some (x,y,z) coordinates from a dataset is:\n" "\n" " 3dmaskdump -mask something+tlrc -o xyzfilename -noijk\n" " '3dcalc( -a dset+tlrc -expr x -datum float )'\n" " '3dcalc( -a dset+tlrc -expr y -datum float )'\n" " '3dcalc( -a dset+tlrc -expr z -datum float )'\n" "\n" " (All of this should be on one command line.)\n" "============================================================================\n" "\n" "Example: 3drotate -prefix Elvis -bshift 10S 0 0 -rotate 30R 0 0 Sinatra+orig\n" "\n" "This will shift the input 10 mm in the superior direction, followed by a 30\n" "degree rotation about the Right-to-Left axis (i.e., nod the head forward).\n" "\n" "============================================================================\n" "Algorithm: The rotation+shift is decomposed into 4 1D shearing operations\n" " (a 3D generalization of Paeth's algorithm). The interpolation\n" " (i.e., resampling) method used for these shears can be controlled\n" " by the following options:\n" "\n" " -Fourier = Use a Fourier method (the default: most accurate; slowest).\n" " -NN = Use the nearest neighbor method.\n" " -linear = Use linear (1st order polynomial) interpolation (least accurate).\n" " -cubic = Use the cubic (3rd order) Lagrange polynomial method.\n" " -quintic = Use the quintic (5th order) Lagrange polynomial method.\n" " -heptic = Use the heptic (7th order) Lagrange polynomial method.\n" "\n" " -Fourier_nopad = Use the Fourier method WITHOUT padding\n" " * If you don't mind - or even want - the wraparound effect\n" " * Works best if dataset grid size is a power of 2, possibly\n" " times powers of 3 and 5, in all directions being altered.\n" " * The main use would seem to be to un-wraparound poorly\n" " reconstructed images, by using a shift; for example:\n" " 3drotate -ashift 30A 0 0 -Fourier_nopad -prefix Anew A+orig\n" " * This option is also available in the Nudge Dataset plugin.\n" "\n" " -clipit = Clip results to input brick range [now the default].\n" " -noclip = Don't clip results to input brick range.\n" "\n" " -zpad n = Zeropad around the edges by 'n' voxels during rotations\n" " (these edge values will be stripped off in the output)\n" " N.B.: Unlike to3d, in this program '-zpad' adds zeros in\n" " all directions.\n" " N.B.: The environment variable AFNI_ROTA_ZPAD can be used\n" " to set a nonzero default value for this parameter.\n" ) ; printf("\n" MASTER_SHORTHELP_STRING ) ; PRINT_COMPILE_DATE ; exit(0) ; } mainENTRY("3drotate main"); machdep(); PRINT_VERSION("3drotate"); AUTHOR("RW Cox"); /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/ { int new_argc ; char ** new_argv ; addto_args( argc , argv , &new_argc , &new_argv ) ; if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; } } AFNI_logger("3drotate",argc,argv) ; #define ERREX(str) (fprintf(stderr,"*** %s\n",str),exit(1)) iopt = 1 ; while( iopt < argc && argv[iopt][0] == '-' ){ if( strncmp(argv[iopt],"-zpad",5) == 0 ){ /* 05 Feb 2001 */ if( zpad > 0 ) fprintf(stderr,"+++ WARNING: second -zpad option!\n") ; zpad = (int) strtod( argv[++iopt] , NULL ) ; if( zpad < 0 ){ fprintf(stderr,"*** ERROR: Can't use -zpad %d\n",zpad) ; exit(1) ; } THD_rota_setpad(zpad,zpad,zpad) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-points",7) == 0 ){ /* 21 Nov 2000 */ dopoints = 1 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-origin",7) == 0 ){ /* 21 Nov 2000 */ xo = strtod( argv[++iopt] , NULL ) ; yo = strtod( argv[++iopt] , NULL ) ; zo = strtod( argv[++iopt] , NULL ) ; doorigin = 1 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-rotpar",7) == 0 ){ /* 07 Feb 2001 */ ATR_float *atr ; if( rotpar_dset != NULL ) ERREX("*** Can't use 2 -rotparent options!") ; if( matvec ) ERREX("*** Can't combine -rotparent with -matvec!") ; if( dcode > 0 || rotarg > 0 ) ERREX("*** Can't use -rotparent with -shift or -rotate options!") ; if( dname != NULL ) ERREX("*** Can't use -rotparent with -dfile or -1Dfile options!") ; rotpar_dset = THD_open_one_dataset( argv[++iopt] ) ; if( rotpar_dset == NULL ) ERREX("*** Can't open -rotparent dataset!\n") ; THD_make_cardinal(rotpar_dset); /* deoblique 21 Oct, 2011 [rickr] */ atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ; if( atr == NULL || atr->nfl < 12 ) ERREX("*** -rotparent dataset doesn't have VOLREG attributes!?") ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-gridpar",7) == 0 ){ /* 07 Feb 2001 */ if( gridpar_dset != NULL ) ERREX("*** Can't use -2 -gridparent options!") ; gridpar_dset = THD_open_one_dataset( argv[++iopt] ) ; if( gridpar_dset == NULL ) ERREX("*** Can't open -gridparent dataset!\n") ; THD_make_cardinal(gridpar_dset); /* deoblique 21 Oct, 2011 [rickr] */ iopt++ ; continue ; } if( strncmp(argv[iopt],"-matvec_",8) == 0 ){ /* 19 Jul 2000 */ MRI_IMAGE *matim ; float *matar , sum ; if( matvec ) ERREX("*** Can't use 2 -matvec options!") ; if( dcode > 0 || rotarg > 0 ) ERREX("*** Can't use -matvec with -shift or -rotate options!") ; if( rotpar_dset != NULL ) ERREX("*** Can't use -matvec with -rotparent option!") ; if( dname != NULL ) ERREX("*** Can't use -matvec with -dfile or -1Dfile options!") ; if( strcmp(argv[iopt],"-matvec_order") == 0 ) matvec = MATVEC_ORDER ; else matvec = MATVEC_DICOM ; if( strcmp(argv[iopt],"-matvec_dset") == 0){ /* 20 July 2000 */ THD_3dim_dataset * mvset ; ATR_float * atr ; mvset = THD_open_dataset( argv[++iopt] ) ; if( mvset == NULL ) ERREX("*** Can't read -matvec_dset dataset!") ; THD_make_cardinal(mvset); /* deoblique 21 Oct, 2011 [rickr] */ atr = THD_find_float_atr( mvset->dblk , "TAGALIGN_MATVEC" ) ; if( atr == NULL || atr->nfl < 12 ) ERREX("*** -matvec_dset doesn't have matrix+vector in .HEAD!") ; matar = atr->fl ; LOAD_DMAT(rmat,matar[0],matar[1],matar[2], matar[4],matar[5],matar[6], matar[8],matar[9],matar[10] ) ; LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ; DSET_delete(mvset) ; } else { matim = mri_read_ascii( argv[++iopt] ) ; if( matim == NULL ) ERREX("Can't read -matvec file!") ; if( matim->nx != 4 || matim->ny != 3 ) ERREX("-matvec file not 4x3!") ; matar = MRI_FLOAT_PTR(matim) ; LOAD_DMAT(rmat,matar[0],matar[1],matar[2], matar[4],matar[5],matar[6], matar[8],matar[9],matar[10] ) ; LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ; mri_free(matim) ; } /* check if matrix is approximately orthogonal */ /* [will be orthogonalized in rot_to_shear_matvec() in thd_shear3d.c] */ pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ; sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0]) +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1]) +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0); if( sum > 0.01 ) ERREX("-matvec matrix not orthogonal!") ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-clipit",4) == 0 ){ /* 11 Apr 2000 */ fprintf(stderr,"++ Notice: -clipit is now the default\n") ; clipit = 1 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-noclip",4) == 0 ){ /* 16 Apr 2002 */ clipit = 0 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-prefix",4) == 0 ){ new_prefix = argv[++iopt] ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-verbose",5) == 0 ){ verb = 1 ; iopt++ ; continue ; } if( strcmp(argv[iopt],"-Fourier_nopad") == 0 ){ /* 13 May 2003 */ THD_rota_method( MRI_FOURIER_NOPAD ) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-Fourier",4) == 0 || strncmp(argv[iopt],"-fourier",4) == 0 ){ THD_rota_method( MRI_FOURIER ) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-cubic",4) == 0 || strncmp(argv[iopt],"-Cubic",4) == 0 ){ THD_rota_method( MRI_CUBIC ) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-quintic",4) == 0 || strncmp(argv[iopt],"-Quintic",4) == 0 ){ THD_rota_method( MRI_QUINTIC ) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-heptic",4) == 0 || strncmp(argv[iopt],"-Heptic",4) == 0 ){ THD_rota_method( MRI_HEPTIC ) ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-linear",4) == 0 || strncmp(argv[iopt],"-Linear",4) == 0 ){ THD_rota_method( MRI_LINEAR ) ; clipit = 0 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-nn",3) == 0 || strncmp(argv[iopt],"-NN",4) == 0 ){ THD_rota_method( MRI_NN ) ; clipit = 0 ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-ashift",4) == 0 ){ if( matvec ) ERREX("*** Can't use -ashift with -matvec!") ; if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);} if( rotpar_dset != NULL ) ERREX("*** Can't use -ashift with -rotparent!") ; if( dname != NULL ) ERREX("*** Can't use -ashift with -dfile or -1Dfile options!") ; dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ; dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ; dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ; dcode = DELTA_AFTER ; iopt++ ; continue ; } if( strncmp(argv[iopt],"-bshift",4) == 0 ){ if( matvec ) ERREX("*** Can't use -bshift with -matvec!") ; if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);} if( rotpar_dset != NULL ) ERREX("*** Can't use -bshift with -rotparent!") ; if( dname != NULL ) ERREX("*** Can't use -bshift with -dfile or -1Dfile options!") ; dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ; dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ; dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ; dcode = DELTA_BEFORE ; iopt++ ; continue ; } #if 0 if( strncmp(argv[iopt],"-cshift",4) == 0 ){ if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);} dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ; dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ; dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ; dcode = DELTA_FIXED ; iopt++ ; continue ; } #endif if( strncmp(argv[iopt],"-rotate",4) == 0 ){ if( matvec ) ERREX("*** Can't use -rotate with -matvec!") ; if( rotarg > 0 ) ERREX("*** Can't have 2 -rotate options!\n") ; if( rotpar_dset != NULL ) ERREX("*** Can't use -rotate with -rotparent!") ; if( dname != NULL ) ERREX("*** Can't use -rotate with -dfile or -1Dfile options!") ; rotarg = iopt ; /* save and process later */ iopt += 4 ; continue ; } /*-- 19 Jun 2001 --*/ if( strcmp(argv[iopt],"-dfile")==0 || strcmp(argv[iopt],"-1Dfile")==0 ){ if( dname != NULL ) ERREX("*** Can't use 2 -dfile/-1Dfile options!") ; if( matvec ) ERREX("*** Can't combine -dfile/-1Dfile with -matvec!") ; if( dcode > 0 || rotarg > 0 ) ERREX("*** Can't use -dfile/-1Dfile with -shift or -rotate options!") ; if( rotpar_dset != NULL ) ERREX("*** Can't use -dfile/-1Dfile with -rotparent!") ; if( strcmp(argv[iopt],"-dfile") ==0 ) dmode = MODE_DFILE ; else if( strcmp(argv[iopt],"-1Dfile")==0 ) dmode = MODE_1DFILE ; dname = argv[++iopt] ; dim = get_dfile_params( dname , dmode ) ; if( dim == NULL ) ERREX("*** Error with -dfile or -1Dfile!") ; dar = MRI_FLOAT_PTR(dim) ; /* 6 x ndar array */ ndar = dim->ny ; iopt++ ; continue ; } fprintf(stderr,"*** Unknown option: %s\n",argv[iopt]) ; exit(1) ; } /*-- check for legal combinations --*/ if( gridpar_dset != NULL && rotpar_dset == NULL ){ fprintf(stderr,"+++ WARNING: -gridparent means nothing without -rotparent!\n"); DSET_delete( gridpar_dset ) ; gridpar_dset = NULL ; } if( gridpar_dset != NULL && dopoints ){ fprintf(stderr,"+++ WARNING: -gridparent means nothing with -points!\n") ; DSET_delete( gridpar_dset ) ; gridpar_dset = NULL ; } if( doorigin && !dopoints ){ fprintf(stderr,"+++ WARNING: -origin means nothing without -points!\n") ; } if( dar != NULL && dopoints ) ERREX("*** Can't combine -dfile/-1Dfile with -points!") ; if( matvec==0 && dcode<0 && rotarg<0 && rotpar_dset==NULL && dar==NULL ) ERREX("Don't you want to do anything [no -rotate,-shift,-matvec,-rotparent,-dfile]?"); /** read input dataset */ if( iopt >= argc && !dopoints ) ERREX("*** No input dataset?") ; if( iopt < argc ){ dset = THD_open_dataset( argv[iopt] ) ; if( dset == NULL ){ fprintf(stderr,"*** Cannot open dataset %s!\n",argv[iopt]); exit(1); } THD_make_cardinal(dset); /* deoblique 21 Oct, 2011 [rickr] */ } else { dset = EDIT_empty_copy(NULL) ; /* 21 Nov 2000: need a fake dataset */ if( !doorigin ) fprintf(stderr,"+++ WARNING: no -origin and no input dataset on command line!\n") ; } if( dopoints && !doorigin ){ xo = dset->daxes->xxorg + 0.5*(dset->daxes->nxx - 1)*dset->daxes->xxdel ; yo = dset->daxes->yyorg + 0.5*(dset->daxes->nyy - 1)*dset->daxes->yydel ; zo = dset->daxes->zzorg + 0.5*(dset->daxes->nzz - 1)*dset->daxes->zzdel ; fprintf(stderr,"+++ Using -origin %g %g %g\n",xo,yo,zo) ; } /* now can process rotation arguments */ ihand = THD_handedness(dset) ; if( rotarg > 0 ){ int neg ; iopt = rotarg ; th1 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ; switch( *cpt ){ default: ERROR_exit("Illegal code after th1 in -rotate\n"); case '\0': case 'x': case 'X': ax1 = 0 ; neg = 0 ; break ; case 'y': case 'Y': ax1 = 1 ; neg = 0 ; break ; case 'z': case 'Z': ax1 = 2 ; neg = 0 ; break ; case 'A': case 'P': case 'R': case 'L': case 'I': case 'S': ax1 = THD_axcode(dset,*cpt) ; neg = (ax1 < 0) ; ax1 = abs(ax1) - 1 ; break ; } if( neg ) th1 = -th1 ; th2 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ; switch( *cpt ){ default: ERROR_exit("Illegal code after th2 in -rotate\n"); case 'x': case 'X': ax2 = 0 ; break ; case '\0': case 'y': case 'Y': ax2 = 1 ; break ; case 'z': case 'Z': ax2 = 2 ; break ; case 'A': case 'P': case 'R': case 'L': case 'I': case 'S': ax2 = THD_axcode(dset,*cpt) ; neg = (ax2 < 0) ; ax2 = abs(ax2) - 1 ; break ; } if( neg ) th2 = -th2 ; th3 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ; switch( *cpt ){ default: ERROR_exit("Illegal code after th3 in -rotate\n"); case 'x': case 'X': ax3 = 0 ; break ; case 'y': case 'Y': ax3 = 1 ; break ; case '\0': case 'z': case 'Z': ax3 = 2 ; break ; case 'A': case 'P': case 'R': case 'L': case 'I': case 'S': ax3 = THD_axcode(dset,*cpt) ; neg = (ax3 < 0) ; ax3 = abs(ax3) - 1 ; break ; } if( neg ) th3 = -th3 ; if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ) INFO_message("All angles after -rotate are 0!") ; if( ax1 < 0 || ax1 > 2 || ax2 < 0 || ax2 > 2 || ax3 < 0 || ax3 > 2 ) ERROR_exit("Can't understand axes codes in -rotate!") ; if( ihand < 0 ){ th1 = -th1 ; th2 = -th2 ; th3 = -th3 ; } #if 0 fprintf(stderr,"ihand=%d th1=%g th2=%g th3=%g\n",ihand,th1,th2,th3); fprintf(stderr,"ax1=%d ax2=%d ax3=%d\n",ax1,ax2,ax3) ; #endif /* notice the minus signs on the angles ------+ */ /* | | | */ if( dopoints ) /* V V V */ rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ; /* 21 Nov 2000 */ } /* may need to process shift arguments as well */ if( dcode > 0 && (cdx != '\0' || cdy != '\0' || cdz != '\0') ){ float qdx=0 , qdy=0 , qdz=0 ; adx = THD_axcode(dset,cdx) ; if( adx > -99 || dx != 0.0 ){ /* 29 Jan 2001: skip if purely 0 */ switch( adx ){ case 1: qdx = -dx ; break ; default: case -1: qdx = dx ; break ; case 2: qdy = -dx ; break ; case -2: qdy = dx ; break ; case 3: qdz = -dx ; break ; case -3: qdz = dx ; break ; } } ady = THD_axcode(dset,cdy) ; if( ady > -99 || dy != 0.0 ){ /* 29 Jan 2001 */ switch( ady ){ case 1: qdx = -dy ; break ; case -1: qdx = dy ; break ; case 2: qdy = -dy ; break ; default: case -2: qdy = dy ; break ; case 3: qdz = -dy ; break ; case -3: qdz = dy ; break ; } } adz = THD_axcode(dset,cdz) ; if( adz > -99 || dz != 0.0 ){ /* 29 Jan 2001 */ switch( adz ){ case 1: qdx = -dz ; break ; case -1: qdx = dz ; break ; case 2: qdy = -dz ; break ; case -2: qdy = dz ; break ; case 3: qdz = -dz ; break ; default: case -3: qdz = dz ; break ; } } if( verb ) fprintf(stderr,"++ Shifting parameters:\n" "++ direction codes: adx=%d ady=%d adz=%d\n" "++ input values: dx=%g dy=%g dz=%g\n" "++ output values: qdx=%g qdy=%g qdz=%g\n" , adx,ady,adz , dx,dy,dz , qdx,qdy,qdz ) ; dx = qdx ; dy = qdy ; dz = qdz ; } /*- 19 July 2000: now can deal with -matvec_dicom case -*/ if( matvec == MATVEC_DICOM ){ /* convert matrix/vector */ pp = DBLE_mat_to_dicomm( dset ) ; /* to dataset coord order */ ppt = TRANSPOSE_DMAT(pp); /* from the DICOM order! */ rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp); tvec = DMATVEC(ppt,tvec); } /*-- 07 Feb 2001: deal with -rotparent and -gridparent case [very similar to -matvec_dicom, actually] --*/ if( gridpar_dset != NULL ){ int mm , nz_gp , nz_ds ; /* check for compatibility! */ mm = THD_dataset_mismatch( gridpar_dset , dset ) ; if( mm & (MISMATCH_DELTA | MISMATCH_ORIENT) ){ fprintf(stderr,"*** Fatal Error:\n" "*** -gridparent dataset and input dataset don't\n" "*** match in grid spacing and/or orientation!\n" ) ; exit(1) ; } if( DSET_NX(gridpar_dset) != DSET_NX(dset) || DSET_NY(gridpar_dset) != DSET_NY(dset) ){ fprintf(stderr,"*** Fatal Error:\n" "*** -gridparent and input datasets\n" "*** don't match in x,y dimensions!\n" ) ; exit(1) ; } /* check for zero padding */ nz_gp = DSET_NZ(gridpar_dset) ; nz_ds = DSET_NZ(dset) ; if( nz_gp < nz_ds ){ fprintf(stderr,"*** Fatal Error:\n" "*** -gridparent has fewer slices than input dataset!\n") ; exit(1) ; } if( nz_gp > nz_ds ){ int npad1 = (nz_gp - nz_ds) / 2 ; int npad2 = (nz_gp - nz_ds) - npad1 ; int add_I=0, add_S=0, add_A=0, add_P=0, add_L=0, add_R=0 ; THD_3dim_dataset * pset ; char *sp1=NULL,*sp2=NULL ; /* where to add slices? and how many? */ switch( dset->daxes->zzorient ){ case ORI_R2L_TYPE: case ORI_L2R_TYPE: add_R=npad1; add_L=npad2; sp1="R"; sp2="L"; break; case ORI_P2A_TYPE: case ORI_A2P_TYPE: add_A=npad1; add_P=npad2; sp1="A"; sp2="P"; break; case ORI_I2S_TYPE: case ORI_S2I_TYPE: add_I=npad1; add_S=npad2; sp1="I"; sp2="S"; break; } /* add them on */ if( verb ) fprintf(stderr,"+++ Zero padding to match -gridparent: -%s %d -%s %d\n", sp1,npad1,sp2,npad2 ) ; pset = THD_zeropad( dset, add_I,add_S,add_A,add_P,add_L,add_R, "Elvis" , ZPAD_PURGE ) ; if( pset == NULL ){ fprintf(stderr,"*** Fatal Error:\n" "*** Can't properly zeropad input dataset!\n" ) ; exit(1) ; } /* toss input datset, replace with padded one */ DSET_delete(dset) ; dset = pset ; } } if( rotpar_dset != NULL ){ /* compute transformation from -rotparent */ ATR_float * atr ; float * matar , sum ; THD_fvec3 fv ; THD_dfvec3 dv,ev,qv , cv_e2, cv_e1, cv_s1, cv_s2 ; /* load transformation from rotparent */ atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ; matar = atr->fl ; LOAD_DMAT(rmat,matar[0],matar[1],matar[2], matar[4],matar[5],matar[6], matar[8],matar[9],matar[10] ) ; LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ; /* check if matrix is orthogonal */ pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ; sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0]) +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1]) +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0); if( sum > 0.01 ) ERREX("-rotparent matrix not orthogonal!") ; /* must alter shift [tvec] to allow for differing coordinates in the rotparent, gridparent, and input datasets */ /* cv_e2 = center of input dataset [Dicom coordinates] */ fv = THD_dataset_center( dset ) ; /* dataset coords */ FVEC3_TO_DFVEC3( fv , cv_e2 ) ; /* convert to double */ /* cv_e1 = center of gridparent */ if( gridpar_dset != NULL ){ fv = THD_dataset_center( gridpar_dset ) ; FVEC3_TO_DFVEC3( fv , cv_e1 ) ; } else { cv_e1 = cv_e2 ; /* what else to do? */ } /* cv_s2 = center of rotation in rotparent */ atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_OLD" ) ; LOAD_DFVEC3( cv_s2 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ; /* cv_s1 = center of base dataset for rotparent */ atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_BASE" ) ; LOAD_DFVEC3( cv_s1 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ; /* compute extra shift due to difference in center of rotation between rotparent and input dataset, then add in shifts caused by -twodup for rotparent and input */ dv = SUB_DFVEC3( cv_e2 , cv_s2 ) ; ev = DMATVEC( rmat , dv ) ; /* R[E2-S2] */ dv = ev ; /* vestige of a stupid bug, since fixed */ ev = SUB_DFVEC3( cv_e1 , cv_s1 ) ; /* E1-S1 */ qv = SUB_DFVEC3( dv , ev ) ; /* R[E2-S2] + S1-E1 */ tvec = ADD_DFVEC3( tvec , qv ) ; /* shifted translation vector */ /* convert transformation from Dicom to dataset coords */ pp = DBLE_mat_to_dicomm( dset ) ; ppt = TRANSPOSE_DMAT(pp); rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp); tvec = DMATVEC(ppt,tvec); /* modify origin of output dataset to match -gridparent */ if( gridpar_dset != NULL ){ dset->daxes->xxorg = gridpar_dset->daxes->xxorg ; dset->daxes->yyorg = gridpar_dset->daxes->yyorg ; dset->daxes->zzorg = gridpar_dset->daxes->zzorg ; /* 12 Feb 2001: adjust origin of time-offsets as well */ if( dset->taxis != NULL && dset->taxis->nsl > 0 ){ dset->taxis->zorg_sl = dset->daxes->zzorg ; } } matvec = MATVEC_ORDER ; /* flag that transform comes from rmat/tvec */ } /*-- 21 Nov 2000: read (x,y,z) points from stdin, process them, quit --*/ if( dopoints ){ THD_dfvec3 xyzorg ; if( !matvec ) LOAD_DFVEC3( tvec , dx,dy,dz ) ; if( dcode < 0 ) dcode = DELTA_AFTER ; LOAD_DFVEC3(xyzorg,xo,yo,zo) ; rotate_stdin_points( xyzorg , rmat , dcode,tvec ) ; /* at end of file */ exit(0) ; } /*-- 12 Feb 2001: adjust time-offsets for slice direction shifts --*/ if( dset->taxis != NULL && dset->taxis->nsl > 0 ){ int ndz ; int kk,jj , nsl = dset->taxis->nsl ; if( matvec ) ndz = (int) rint( tvec.xyz[2] / fabs(dset->daxes->zzdel) ) ; /* shift */ else ndz = (int) rint( dz / fabs(dset->daxes->xxdel) ) ; if( ndz != 0 ){ float * tsl = (float *)malloc(sizeof(float)*nsl) ; for( kk=0 ; kk < nsl ; kk ++ ){ jj = kk - ndz ; if( jj < 0 || jj >= nsl ) tsl[kk] = 0.0 ; else tsl[kk] = dset->taxis->toff_sl[jj] ; } EDIT_dset_items( dset , ADN_toff_sl , tsl , ADN_none ) ; free(tsl) ; if( verb ) fprintf(stderr,"+++ adjusting time-offsets by %d slices\n",ndz) ; } } /*- read dataset, prepare to process it, write back out (with new name) -*/ DSET_mallocize(dset) ; DSET_load(dset) ; dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */ dset->idcode = MCW_new_idcode() ; /* 08 Jun 1999 - is a new dataset */ EDIT_dset_items( dset , ADN_prefix , new_prefix , ADN_label1 , new_prefix , ADN_none ) ; if( THD_deathcon() && THD_is_file(dset->dblk->diskptr->header_name) ){ fprintf(stderr, "** ERROR: Output file %s already exists -- cannot continue!\n", dset->dblk->diskptr->header_name ) ; exit(1) ; } /* old history is already in the dataset */ tross_Make_History( "3drotate" , argc,argv , dset ) ; nvox = DSET_NVOX(dset) ; fvol = (float *) malloc( sizeof(float) * nvox ) ; nval = DSET_NVALS(dset) ; if( verb ){ fprintf(stderr,"+++ %d sub-bricks: ",nval) ; cputim = COX_cpu_time() ; } /* 03 May 2005: save rotation center */ { THD_fvec3 cv ; cv = THD_dataset_center( dset ) ; THD_set_float_atr( dset->dblk , "ROTATE_CENTER_OLD" , 3 , cv.xyz ) ; THD_set_float_atr( dset->dblk , "ROTATE_CENTER_BASE" , 3 , cv.xyz ) ; } /*-- loop over all sub-bricks: copy into fvol, rotate fvol, copy back --*/ for( ival=0 ; ival < nval ; ival++ ){ if( verb ) fprintf(stderr,"%d",ival) ; EDIT_coerce_type( nvox , DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) , MRI_float,fvol ) ; if( verb ) fprintf(stderr,".") ; if( clipit ){ /* 11 Apr 2000 */ register int ii ; register float bb,tt ; bb = tt = fvol[0] ; for( ii=1 ; ii < nvox ; ii++ ){ if( fvol[ii] < bb ) bb = fvol[ii] ; else if( fvol[ii] > tt ) tt = fvol[ii] ; } cbot = bb ; ctop = tt ; } /* 19 Jun 2001: get matrix/vector from rot params in file */ skipit = 0 ; if( dar != NULL ){ THD_dvecmat dvm ; char rotcom[256] ; int jj=ival ; static int ndar_over=0 ; if( jj >= ndar ){ jj = ndar-1 ; if( ndar_over == 0 ) WARNING_message("from brick %d on, using last line (%d) from %s\n", jj , ndar-1 , dname ) ; ndar_over++ ; } sprintf(rotcom,"-rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP", dar[0+6*jj] , dar[1+6*jj] , dar[2+6*jj] , dar[3+6*jj] , dar[4+6*jj] , dar[5+6*jj] ) ; dvm = THD_rotcom_to_matvec( dset , rotcom ) ; /* thd_rotangles.c */ rmat = dvm.mm ; tvec = dvm.vv ; matvec = MATVEC_ORDER ; skipit = (dar[0+6*jj]==0.0 && dar[1+6*jj]==0.0 && dar[2+6*jj]==0.0 && dar[3+6*jj]==0.0 && dar[4+6*jj]==0.0 && dar[5+6*jj]==0.0 ); if( !skipit ){ skipit = ( fabs(rmat.mat[0][0]-1.0) < 0.00001 ) && ( fabs(rmat.mat[1][1]-1.0) < 0.00001 ) && ( fabs(rmat.mat[2][2]-1.0) < 0.00001 ) && ( fabs(tvec.xyz[0]) < 0.001 ) && ( fabs(tvec.xyz[1]) < 0.001 ) && ( fabs(tvec.xyz[2]) < 0.001 ) ; } if( skipit ) fprintf(stderr,"[Matrix near identity, skipping]"); } /** carry out the rotation **/ if( !skipit ){ if( matvec ){ THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) , fabs(DSET_DX(dset)) , fabs(DSET_DY(dset)) , fabs(DSET_DZ(dset)) , fvol , rmat , tvec ) ; } else { THD_rota_vol( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) , fabs(DSET_DX(dset)), fabs(DSET_DY(dset)), fabs(DSET_DZ(dset)), fvol , ax1,th1, ax2,th2, ax3,th3, dcode,dx,dy,dz ) ; /* 02 May 2005: save matrix+vector for recording below */ rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ; LOAD_DFVEC3( tvec , dx,dy,dz ) ; if( dcode == DELTA_BEFORE ) tvec = DMATVEC(rmat,tvec) ; } } /** 02 May 2005: record operation in header **/ { float matar[12] ; THD_dmat33 drmat ; THD_dfvec3 dvec ; char anam[64] ; /* must transform [rmat,tvec] from dataset to DICOM coords */ /* (the inverse of the DICOM to dataset transforms earlier) */ pp = DBLE_mat_to_dicomm(dset) ; ppt = TRANSPOSE_DMAT(pp) ; drmat = DMAT_MUL(pp,rmat) ; drmat = DMAT_MUL(drmat,ppt); dvec = DMATVEC (pp,tvec) ; UNLOAD_DMAT(drmat,matar[0],matar[1],matar[2], matar[4],matar[5],matar[6], matar[8],matar[9],matar[10] ) ; UNLOAD_DFVEC3(dvec,matar[3],matar[7],matar[11]) ; sprintf(anam,"ROTATE_MATVEC_%06d",ival) ; THD_set_float_atr( dset->dblk , anam , 12 , matar ) ; } if( verb ) fprintf(stderr,".") ; if( clipit ){ /* 11 Apr 2000 */ register int ii ; register float bb,tt ; bb = cbot ; tt = ctop ; for( ii=0 ; ii < nvox ; ii++ ){ if( fvol[ii] < bb ) fvol[ii] = bb ; else if( fvol[ii] > tt ) fvol[ii] = tt ; } } EDIT_coerce_type( nvox , MRI_float,fvol , DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ) ; } /* end of loop over sub-bricks */ if( verb ){ cputim = COX_cpu_time() - cputim ; fprintf(stderr,"\n+++ CPU time=%10.3g s" , cputim) ; if( nval > 1 ) fprintf(stderr," [= %10.3g s/sub-brick]" , cputim/nval) ; fprintf(stderr,"\n+++ Writing dataset to disk in %s", DSET_BRIKNAME(dset) ) ; } dset->dblk->master_nvals = 0 ; /* 11 Apr 2000 hack */ DSET_write(dset) ; if( verb ) fprintf(stderr,"\n") ; exit(0) ; } /*=================================================================================*/ #include #define NBUF 1024 /* line buffer size */ #define DUPOUT(n) fprintf(fpout,"%s",linbuf+n) static void rotate_stdin_points( THD_dfvec3 xyzorg, THD_dmat33 rmat, int dcode, THD_dfvec3 tvec ) { char linbuf[NBUF] , *cp ; FILE *fpin=stdin , *fpout=stdout ; int ii , kk , nbuf , ll , nn , ld ; double xx,yy,zz ; THD_dfvec3 xyz ; if( verb ){ DUMP_DMAT33("Rotation",rmat) ; } /** rmat = DMAT_INV(rmat) ; **/ /*-- loop over input lines --*/ ll = ld = 0 ; while(1){ ll++ ; /* line count */ cp = afni_fgets( linbuf , NBUF , fpin ) ; /* read the line */ if( cp == NULL ) break ; /* end of file => end of loop */ kk = strlen(linbuf) ; if( kk == 0 ) continue ; /* empty line => get next line */ /* find 1st nonblank */ for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ; /* nada */ if( ii == kk || /* all blanks */ (linbuf[ii] == '/' && linbuf[ii+1] == '/') ){ /* or comment */ DUPOUT(0) ; continue ; } /* scan line for data */ nn = sscanf(linbuf+ii , "%lf%lf%lf%n" , &xx,&yy,&zz,&nbuf ) ; if( nn < 3 ){ fprintf(stderr,"+++ WARNING: input line %d was incomplete\n",ll) ; continue ; } nbuf += ii ; /* position of next character after zz */ /* process vector */ LOAD_DFVEC3(xyz , xx,yy,zz) ; xyz = SUB_DFVEC3(xyz,xyzorg) ; if( dcode == DELTA_BEFORE ) xyz = ADD_DFVEC3(xyz,tvec) ; xyz = DMATVEC(rmat,xyz) ; if( dcode == DELTA_AFTER ) xyz = ADD_DFVEC3(xyz,tvec) ; xyz = ADD_DFVEC3(xyz,xyzorg) ; fprintf(fpout,"%g %g %g%s",xyz.xyz[0],xyz.xyz[1],xyz.xyz[2],linbuf+nbuf) ; ld++ ; } /* end of loop over input lines */ if( verb ) fprintf(stderr,"-points: read %d lines, wrote %d lines\n",ll-1,ld) ; } /*----------------------------------------------------------------------- 19 Jun 2001: Get a 6 x N image of the rotation parameters -------------------------------------------------------------------------*/ MRI_IMAGE * get_dfile_params( char * fname , int mode ) { MRI_IMAGE *outim , *flim ; float *oar , *far ; int nx,ny , ii,jj ; if( fname == NULL || fname[0] == '\0' ) return NULL ; flim = mri_read_ascii( fname ) ; if( flim == NULL ) return NULL ; nx = flim->nx ; ny = flim->ny ; far = MRI_FLOAT_PTR(flim) ; if( mode == MODE_DFILE ){ /* skip 1st element of each row */ if( nx < 7 ){ fprintf(stderr,"** -dfile %s has too few columns!\n",fname) ; mri_free(flim) ; return NULL ; } outim = mri_new( 6 , ny , MRI_float ) ; oar = MRI_FLOAT_PTR(outim) ; for( jj=0 ; jj < ny ; jj++ ){ oar[0+6*jj] = far[1+nx*jj] ; oar[1+6*jj] = far[2+nx*jj] ; oar[2+6*jj] = far[3+nx*jj] ; oar[3+6*jj] = far[4+nx*jj] ; oar[4+6*jj] = far[5+nx*jj] ; oar[5+6*jj] = far[6+nx*jj] ; } mri_free(flim) ; } else if( mode == MODE_1DFILE ){ /* first 6 elements of each row */ if( nx < 6 ){ fprintf(stderr,"** -1Dfile %s has too few columns!\n",fname) ; mri_free(flim) ; return NULL ; } else if( nx == 6 ){ outim = flim ; } else { outim = mri_new( 6 , ny , MRI_float ) ; oar = MRI_FLOAT_PTR(outim) ; for( jj=0 ; jj < ny ; jj++ ){ oar[0+6*jj] = far[0+nx*jj] ; oar[1+6*jj] = far[1+nx*jj] ; oar[2+6*jj] = far[2+nx*jj] ; oar[3+6*jj] = far[3+nx*jj] ; oar[4+6*jj] = far[4+nx*jj] ; oar[5+6*jj] = far[5+nx*jj] ; } mri_free(flim) ; } } else { fprintf(stderr,"** get_dfile_params: illegal mode=%d\n",mode) ; mri_free(flim) ; return NULL ; } return outim ; }