/***************************************************************************** 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 "mrilib.h" /*-- these macros stolen from file thd.h --*/ #define ORCODE(aa) \ ( (aa)=='R' ? ORI_R2L_TYPE : (aa)=='L' ? ORI_L2R_TYPE : \ (aa)=='P' ? ORI_P2A_TYPE : (aa)=='A' ? ORI_A2P_TYPE : \ (aa)=='I' ? ORI_I2S_TYPE : (aa)=='S' ? ORI_S2I_TYPE : ILLEGAL_TYPE ) #define OR3OK(x,y,z) ( ((x)&6) + ((y)&6) + ((z)&6) == 6 ) /*------------------------------------------------------------------------*/ void Syntax(int detail) { printf( "Usage: 3dUndump [options] infile ...\n" "Assembles a 3D dataset from an ASCII list of coordinates and\n" "(optionally) values.\n" "\n" "Options:\n" "--------\n" " -prefix ppp = 'ppp' is the prefix for the output dataset\n" " [default = undump].\n" "\n" " -master mmm = 'mmm' is the master dataset, whose geometry\n" " *OR* will determine the geometry of the output.\n" " -dimen I J K = Sets the dimensions of the output dataset to\n" " be I by J by K voxels. (Each I, J, and K\n" " must be >= 1.) This option can be used to\n" " create a dataset of a specific size for test\n" " purposes, when no suitable master exists.\n" " ** N.B.: Exactly one of -master or -dimen must be given.\n" " ==>> Unless -ROImask is used!\n" "\n" " -mask MMM = This option specifies a mask dataset 'MMM', which\n" " will control which voxels are allowed to get\n" " values set. If the mask is present, only\n" " voxels that are nonzero in the mask can be\n" " set in the new dataset.\n" " * A mask can be created from a pre-existing dataset\n" " with program 3dAutomask.\n" " * Combining a mask with sphere insertion makes\n" " a lot of sense (to me, at least).\n" "\n" " -datum type = 'type' determines the voxel data type of the\n" " output, which may be byte, short, or float\n" " [default = short].\n" "\n" " -dval vvv = 'vvv' is the default value stored in each\n" " input voxel that does not have a value\n" " supplied in the input file [default = 1].\n" " *** For example: '-dval 7' makes all the voxels\n" " whose locations (via '-ijk' or '-xyz') are\n" " given without a value get the value 7.\n" "\n" " -fval fff = 'fff' is the fill value, used for each voxel\n" " in the output dataset that is NOT listed\n" " in the input file [default = 0].\n" " *** For example: '-fval 1' makes all the completely\n" " unspecified voxels get the value 1.\n" "\n" " -ijk = Coordinates in the input file are (i,j,k) index\n" " *OR* triples, as might be output by 3dmaskdump.\n" " -xyz = Coordinates in the input file are (x,y,z)\n" " spatial coordinates, in mm. If neither\n" " -ijk or -xyz is given, the default is -ijk.\n" " ** N.B.: -xyz can only be used with -master. If -dimen\n" " is used to specify the size of the output dataset,\n" " (x,y,z) coordinates are not defined (until you\n" " use 3drefit to define the spatial structure).\n" " ** N.B.: Also see the -ROImask option (infra) for another\n" " way to specify what voxels in the output dataset\n" " get what values:\n" " -- with -ROImask, neither -xyz nor -ijk is used.\n" "\n" " -srad rrr = Specifies that a sphere of radius 'rrr' will be\n" " filled about each input (x,y,z) or (i,j,k) voxel.\n" " If the radius is not given, or is 0, then each\n" " input data line sets the value in only one voxel.\n" " * If '-master' is used, then 'rrr' is in mm.\n" " * If '-dimen' is used, then 'rrr' is in voxels.\n" " *** For example: '-srad 7' means put a sphere of\n" " radius 7 mm about each input point.\n\n" " -cubes = Put cubes down instead of spheres. The \"radius\" then\n" " is half the length of a side.\n" "\n" " -orient code = Specifies the coordinate order used by -xyz.\n" " The code must be 3 letters, one each from the pairs\n" " {R,L} {A,P} {I,S}. The first letter gives the\n" " orientation of the x-axis, the second the orientation\n" " of the y-axis, the third the z-axis:\n" " R = right-to-left L = left-to-right\n" " A = anterior-to-posterior P = posterior-to-anterior\n" " I = inferior-to-superior S = superior-to-inferior\n" " If -orient isn't used, then the coordinate order of the\n" " -master dataset is used to interpret (x,y,z) inputs.\n" " ** N.B.: If -dimen is used (which implies -ijk), then the\n" " only use of -orient is to specify the axes ordering\n" " of the output dataset. If -master is used instead,\n" " the output dataset's axes ordering is the same as the\n" " -master dataset's, regardless of -orient.\n" " ** You probably don't need this option, and it is kept\n" " here only for historical purposes.\n" "\n" " -head_only = A 'secret' option for creating only the .HEAD file which\n" " gets exploited by the AFNI matlab library function\n" " New_HEAD.m\n" "\n" " -ROImask rrr = This option that lets you specify which voxels get what\n" " numbers by using a dataset 'rrr', instead of coordinates.\n" " ==>>** With this method, the input file should have just\n" " one number per line (trailing numbers will be ignored).\n" " ** Due to the special way that 3dUndump reads input files, you\n" " CANNOT specify an input file using the 1D '[subscript]'\n" " notation to pick out a single column of a multicolumn\n" " file. Instead, you can do something like\n" " 1dcat file.1D'[3]' | 3dUndump -ROImask rmask+orig -prefix ppp -\n" " where the last '-' says to read from standard input.\n" " ** A more complicated example, using an ROI mask dataset 'mmm+orig'\n" " to extract average values from a functional dataset, then create\n" " a dataset where the values stored are the ROI averages:\n" " 3dROIstats -mask mmm+orig -1Dformat func+orig'[1]' | \\\n" " | 1dcat stdin: | 3dUndump -prefix uuu -datum float -ROImask mmm+orig -\n" " Again, the final '-' tells 3dUndump to read the values to be\n" " stored from standard input (the pipe).\n" " ** If the numbers in the input file are fractional (e.g., '1.372'),\n" " be sure to use the '-datum float' option -- otherwise, the\n" " default output is '-datum short', which will truncate values!\n" " * The 'rrr' dataset must be of integer type -- that is,\n" " the values inside must be bytes or shorts. If you don't\n" " know, use program 3dinfo to check.\n" " * All voxels with value 1 in dataset 'rrr' get the number in\n" " the first row of the input file.\n" " * All voxels with value 2 in dataset 'rrr' get the number in\n" " the second row of the input file.\n" " * Et cetera -- all voxels with value 'n' in dataset 'rrr' get\n" " the number in the n-th row of the input file.\n" " * Zero or negative values in 'rrr' are ignored completely.\n" " * The output dataset has the same spatial grid as 'rrr'\n" " (i.e., as if '-master rrr' were used).\n" " * The following options cannot be used with -ROImask:\n" " -dimen -master -mask -dval\n" " -ijk -xyz -srad -orient -head_only\n" " * This option was added 09 Nov 2011:\n" " -- Happy 280th Birthday to Benjamin Banneker!\n" " -- http://en.wikipedia.org/wiki/Benjamin_Banneker\n" "\n" "Input File Format:\n" "------------------\n" " The input file(s) are ASCII files, with one voxel specification per\n" " line. A voxel specification is 3 numbers (-ijk or -xyz coordinates),\n" " with an optional 4th number giving the voxel value. For example:\n" "\n" " 1 2 3 \n" " 3 2 1 5\n" " 5.3 6.2 3.7\n" " // this line illustrates a comment\n" "\n" " The first line puts a voxel (with value given by '-dval') at point\n" " (1,2,3). The second line puts a voxel (with value 5) at point (3,2,1).\n" " The third line puts a voxel (with value given by '-dval') at point\n" " (5.3,6.2,3.7). If -ijk is in effect, and fractional coordinates\n" " are given, they will be rounded to the nearest integers; for example,\n" " the third line would be equivalent to (i,j,k) = (5,6,4).\n" "\n" "Notes:\n" "------\n" "* This program creates a 1 sub-brick file. You can 'glue' multiple\n" " files together using 3dbucket or 3dTcat to make multi-brick datasets.\n" " *** At this time, 3dUndump cannot create a multi-brick dataset :-(\n" "\n" "* If one input filename is '-', then stdin will be used for input.\n" " This feature is for the intrepid Unix user who wants to pipe the\n" " input into 3dUndump from another program.\n" "\n" "* If no input files are given, an 'empty' dataset is created.\n" " For example, to create an all zero 3D dataset with 1 million voxels:\n" " 3dUndump -dimen 100 100 100 -prefix AllZeroAFNI\n" " 3dUndump -dimen 100 100 100 -prefix AllZeroNIFTI.nii\n" " *** This is probably the simplest way to create an all zero dataset\n" " with given grid dimensions in AFNI, without any pre-existing\n" " 'master' dataset to start with. If you want to further change\n" " the voxel sizes (in mm), you can use 3drefit to alter such\n" " parameters after the initial act of creation ab nihilo.\n" " *** You can combine 3dUndump with 3dcalc to create an all zero\n" " 3D+time dataset from 'thin air', as in\n" " 3dUndump -dimen 128 128 32 -prefix AllZero_A -datum float\n" " 3dcalc -a AllZero_A+orig -b '1D: 100@' -expr 0 -prefix AllZero_B\n" " If you replace the '0' expression with 'gran(0,1)', you'd get a\n" " random 3D+time dataset, which might be useful for testing purposes.\n" "\n" "* By default, the output dataset is of type '-fim', unless the -master\n" " dataset is an anat type. You can change the output type later using 3drefit.\n" "\n" "* You could use program 1dcat to extract specific columns from a\n" " multi-column rectangular file (e.g., to get a specific sub-brick\n" " from the output of 3dmaskdump), and use the output of 1dcat as input\n" " to this program. If you know what you are doing, that is.\n" "\n" "* [19 Feb 2004] The -mask and -srad options were added this day.\n" " Also, a fifth value on an input line, if present, is taken as a\n" " sphere radius to be used for that input point only. Thus, input\n" " 3.3 4.4 5.5 6.6 7.7\n" " means to put the value 6.6 into a sphere of radius 7.7 mm centered\n" " about (x,y,z)=(3.3,4.4,5.5).\n" "\n" "* [10 Nov 2008] Commas (',') inside an input line are converted to\n" " spaces (' ') before the line is interpreted. This feature is for\n" " convenience for people writing files in CSV (Comma Separated Values)\n" " format.\n" " ++ [14 Feb 2010] Semicolons (';') and colons (':') are now changed\n" " to blanks, as well. In addition, any line that starts with\n" " an alphabetic character, or with '#' or '/' will be skipped\n" " (presumably it is some kind of comment).\n" "\n" "* [31 Dec 2008] Inputs of 'NaN' are explicitly converted to zero, and\n" " a warning message is printed. AFNI programs do not like with NaN\n" " floating point values!\n" "\n" "* [09 Jun 2021] The new option '-allow_NaN' will let NaN (not-a-number)\n" " values be entered. This option must be used before any input values\n" " are given (e.g., via '-dval').\n" " ++ NaN is a floating point concept. If you input (say) -dval as NaN,\n" " but then use it in creating a short-valued dataset, you will get\n" " the same as '-dval 0' -- there is no equivalent to Not-a-Number\n" " in the 2's complement world of integer representations.\n" " ++ IN OTHER WORDS: use '-datum float' if you use '-allow_Nan',\n" " OR you will end up mystified.\n" " ++ Please note that AFNI programs will convert NaN values to 0 on input,\n" " from files stored in NIfTI and some other more obscure formats.\n" " ++ And: .BRIK files will be similarly scanned/fixed on input if\n" " environment variable AFNI_FLOATSCAN is set to YES.\n" " ++ And: the AFNI GUI will not like you at some point if you are\n" " trying to view a file containing NaNs.\n" " ++ In other words:\n" " The only reason for '-allow_NaN' is to create files for testing.\n" " And only for the semi-deranged amongst us. (You know who I mean.)\n" " ++ Sample command line:\n" " echo '0 0 0' | 3dUndump -allow_NaN -datum float -dimen 50 50 50 -prefix Ubad -dval NaN -\n" " float_scan Ubad+orig.BRIK\n" " Program float_scan should tell you that Ubad+orig.BRIK has 1 bad value in it.\n" "\n" "-- RWCox -- October 2000\n" ) ; PRINT_COMPILE_DATE ; return ; } /*---------------------------------------------------------------------------*/ #define NBUF 1024 /* line buffer size */ int main( int argc , char * argv[] ) { int do_ijk=1 , dimen_ii=0 , dimen_jj=0 , dimen_kk=0 , datum=MRI_short ; THD_3dim_dataset *mset=NULL ; char *prefix="undump" , *orcode=NULL ; THD_coorder cord ; float dval_float=1.0f, fval_float=0.0f, *fbr=NULL ; short dval_short=1 , fval_short=0 , *sbr=NULL ; byte dval_byte =1 , fval_byte =0 , *bbr=NULL , *mmask=NULL ; FILE *fp ; THD_3dim_dataset *dset , *maskset=NULL ; int iarg , ii,jj,kk,ll,ijk , nx,ny,nz , nxyz , nn , do_head_only=0; float xx,yy,zz,vv=0.0f; short sv=0 ; byte bv=0 ; char linbuf[NBUF] , *cp ; byte cubes = 0; float xxdown,xxup , yydown,yyup , zzdown,zzup ; float srad=0.0f, vrad,rii,rjj,rkk,qii,qjj,qkk , dx,dy,dz ; /* 19 Feb 2004 */ THD_3dim_dataset *ROImask=NULL ; /* 09 Nov 2011 */ int have_dimen=0, have_master=0, have_mask=0, have_dval=0, have_ijk=0, have_xyz=0, have_srad=0, have_orient=0, have_head_only=0 ; short ROInum , *ROIsss=NULL ; int nfill=0 ; int allow_NaN=0 ; /* 09 Jun 2021 */ #define MAYBE_SCAN(a,b) ( (allow_NaN) ? (0) : thd_floatscan((a),(b)) ) /*-- help? --*/ /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/ mainENTRY("3dUndump main") ; machdep() ; PRINT_VERSION("3dUndump"); machdep() ; { 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("3dUndump",argc,argv) ; /*-- command line options --*/ do_head_only = 0 ; iarg = 1 ; while( iarg < argc && argv[iarg][0] == '-' ){ if( strcmp(argv[iarg],"-help") == 0 || strcmp(argv[iarg],"-h") == 0 ) { Syntax(strlen(argv[iarg])>3 ? 2:1) ; exit(0); } if( strcmp(argv[iarg],"-") == 0 ){ /* a single - is an input filename */ break ; } if( strcmp(argv[iarg],"-allow_NaN") == 0 ){ /* 09 Jun 2021 */ allow_NaN = 1 ; iarg++ ; continue ; } /*-----*/ if( strcasecmp(argv[iarg],"-ROImask") == 0 ){ /* 09 Nov 2011 */ if( ROImask != NULL ) ERROR_exit("Can't have 2 -ROImask options!") ; if( iarg+1 >= argc ) ERROR_exit("-ROImask: no argument follows!?") ; ROImask = THD_open_dataset( argv[++iarg] ) ; if( ROImask == NULL ) ERROR_exit("-ROImask: can't open dataset '%s'",argv[iarg] ) ; ii = (int)DSET_BRICK_TYPE(ROImask,0) ; if( ii != MRI_byte && ii != MRI_short ) ERROR_exit("-ROImask: dataset values are '%s', not 'byte' or 'short'", MRI_TYPE_name[ii] ) ; DSET_load(ROImask) ; CHECK_LOAD_ERROR(ROImask) ; ROIsss = (short *)malloc( sizeof(short)*DSET_NVOX(ROImask) ) ; EDIT_coerce_type( DSET_NVOX(ROImask) , ii , DSET_BRICK_ARRAY(ROImask,0) , MRI_short , ROIsss ) ; DSET_unload(ROImask) ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-prefix") == 0 ){ if( iarg+1 >= argc ) ERROR_exit("-prefix: no argument follows!?") ; else if( !THD_filename_ok(argv[++iarg]) ) ERROR_exit("-prefix: Illegal prefix given!") ; prefix = argv[iarg] ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-master") == 0 ){ if( iarg+1 >= argc ) ERROR_exit("-master: no argument follows!?") ; else if( mset != NULL ) ERROR_exit("-master: can't have 2 -master options!") ; else if( dimen_ii > 0 ) ERROR_exit("-master: conflicts with previous -dimen!") ; mset = THD_open_dataset( argv[++iarg] ) ; if( mset == NULL ) ERROR_exit("-master: can't open dataset" ) ; have_master++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-mask") == 0 ){ if( iarg+1 >= argc ) ERROR_exit("-mask: no argument follows!?") ; else if( maskset != NULL ) ERROR_exit("-mask: can't have two -mask options!") ; maskset = THD_open_dataset( argv[++iarg] ) ; if( maskset == NULL ) ERROR_exit("-mask: can't open dataset" ) ; have_mask++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-dimen") == 0 ){ if( iarg+3 >= argc ) ERROR_exit("-dimen: don't have 3 arguments following!?") ; else if( mset != NULL ) ERROR_exit("-dimen: conflicts with previous -master!") ; else if( dimen_ii > 0 ) ERROR_exit("-dimen: can't have two -dimen options!") ; dimen_ii = strtol(argv[++iarg],NULL,10) ; dimen_jj = strtol(argv[++iarg],NULL,10) ; dimen_kk = strtol(argv[++iarg],NULL,10) ; if( dimen_ii < 1 || dimen_jj < 1 || dimen_kk < 1 ) ERROR_exit("-dimen: values following are not all >= 1!") ; have_dimen++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-datum") == 0 ){ if( ++iarg >= argc ) ERROR_exit("-datum: no argument follows?!") ; if( strcmp(argv[iarg],"short") == 0 ) datum = MRI_short ; else if( strcmp(argv[iarg],"float") == 0 ) datum = MRI_float ; else if( strcmp(argv[iarg],"byte") == 0 ) datum = MRI_byte ; else ERROR_exit("-datum: illegal type given!") ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-srad") == 0 ){ /* 19 Feb 2004 */ if( iarg+1 >= argc ) ERROR_exit("-srad: no argument follows!?") ; srad = strtod( argv[++iarg] , NULL ) ; if( srad <= 0.0 || thd_floatscan(1,&srad) ){ WARNING_message("-srad '%s' is ignored!",argv[iarg]); srad = 0.0 ; } have_srad++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-cubes") == 0 ){ /* 19 Feb 2004 */ cubes = 1 ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-dval") == 0 ){ if( iarg+1 >= argc ) ERROR_exit("-dval: no argument follows!?") ; dval_float = strtod( argv[++iarg] , NULL ) ; if( MAYBE_SCAN(1,&dval_float) ) ERROR_exit("Illegal value '%s'entered for -dval!",argv[iarg]) ; if( isnan(dval_float) ){ WARNING_message("NaN value being used for -dval") ; dval_short = 0 ; dval_byte = 0 ; } else { dval_short = (short) rintf(dval_float) ; dval_byte = (byte) dval_short ; } have_dval++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-fval") == 0 ){ if( iarg+1 >= argc ) ERROR_exit("-fval: no argument follows!?") ; fval_float = strtod( argv[++iarg] , NULL ) ; if( MAYBE_SCAN(1,&fval_float) ) ERROR_exit("Illegal value '%s'entered for -fval!",argv[iarg]) ; if( isnan(fval_float) ){ WARNING_message("NaN value being used for -fval") ; fval_short = 0 ; fval_byte = 0 ; } else { fval_short = (short) rintf(fval_float) ; fval_byte = (byte) dval_short ; } iarg++ ; continue ; } if( strcmp(argv[iarg],"-ijk") == 0 ){ do_ijk = 1 ; have_ijk++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-xyz") == 0 ){ do_ijk = 0 ; have_xyz++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-head_only") == 0 ){ do_head_only = 1 ; have_head_only++ ; iarg++ ; continue ; } /*-----*/ if( strcmp(argv[iarg],"-orient") == 0 ){ int xx,yy,zz ; if( iarg+1 >= argc ) ERROR_exit("-orient: no argument follows!?") ; orcode = argv[++iarg] ; if( strlen(orcode) != 3 ) ERROR_exit("-orient: illegal argument follows") ; xx = ORCODE(orcode[0]) ; yy = ORCODE(orcode[1]) ; zz = ORCODE(orcode[2]) ; if( xx < 0 || yy < 0 || zz < 0 || !OR3OK(xx,yy,zz) ) ERROR_exit("-orient: illegal argument follows") ; have_orient++ ; iarg++ ; continue ; } /*-----*/ ERROR_message("Unknown option: %s",argv[iarg]) ; suggest_best_prog_option(argv[0], argv[iarg]); exit(1); } /* end of loop over command line options */ if( argc < 3) { ERROR_message("Too few options"); Syntax(0) ; exit(1); } /*-- check for inconsistencies --*/ #if 0 if( iarg >= argc && !do_head_only) ERROR_exit("No input files on command line!?") ; #endif if( datum != MRI_float && allow_NaN ){ /* 10 Jun 2021 */ WARNING_message("'-allow_NaN' without '-datum float' is meaningless!") ; if( isnan(fval_float) ){ fval_byte = 0; fval_short = 0; fval_float = 0.0f; } if( isnan(dval_float) ){ dval_byte = 0; dval_short = 0; dval_float = 0.0f; } } if( ROImask != NULL && ( have_dimen || have_master || have_mask || have_dval || have_ijk || have_xyz || have_srad || have_orient || have_head_only ) ) ERROR_exit("-ROImask conflicts with one (or more) of these options:\n" " -dimen -master -mask -dval -ijk\n" " -xyz -srad -orient -head_only\n" ) ; if( do_ijk == 0 && mset == NULL && ROImask == NULL ) ERROR_exit("Can't use -xyz without -master also!") ; if( mset == NULL && ROImask == NULL && dimen_ii < 2 ) ERROR_exit("Must use exactly one of -master or -dimen options on command line"); if( (datum == MRI_short && dval_short == fval_short) || (datum == MRI_float && dval_float == fval_float) || (datum == MRI_byte && dval_byte == fval_byte ) ){ WARNING_message("-dval and -fval are the same!") ; } /*-- set orcode to value from -master, if this is needed --*/ if( mset != NULL && ROImask == NULL && do_ijk == 0 && orcode == NULL ){ orcode = malloc(4) ; orcode[0] = ORIENT_typestr[mset->daxes->xxorient][0] ; orcode[1] = ORIENT_typestr[mset->daxes->yyorient][0] ; orcode[2] = ORIENT_typestr[mset->daxes->zzorient][0] ; orcode[3] = '\0' ; } THD_coorder_fill( orcode , &cord ) ; /* setup coordinate order */ /*-- make empty dataset --*/ if( mset != NULL || ROImask != NULL ){ /* from -master or -ROImask */ THD_3dim_dataset *qset = (mset != NULL) ? mset : ROImask ; dset = EDIT_empty_copy( qset ) ; EDIT_dset_items( dset , ADN_prefix , prefix , ADN_datum_all , datum , ADN_nvals , 1 , ADN_ntt , 0 , ADN_func_type , ISANAT(qset) ? qset->func_type : FUNC_FIM_TYPE , ADN_directory_name , "./" , ADN_none ) ; } else { /* from nothing */ THD_ivec3 iv_nxyz , iv_xyzorient ; THD_fvec3 fv_xyzorg , fv_xyzdel ; LOAD_IVEC3( iv_nxyz , dimen_ii , dimen_jj , dimen_kk ) ; LOAD_IVEC3( iv_xyzorient , cord.xxor , cord.yyor , cord.zzor ) ; LOAD_FVEC3( fv_xyzdel , ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? 1.0 : -1.0 , ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? 1.0 : -1.0 , ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? 1.0 : -1.0 ) ; LOAD_FVEC3( fv_xyzorg , ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? -0.5*dimen_ii : 0.5*dimen_ii, ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? -0.5*dimen_jj : 0.5*dimen_jj, ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? -0.5*dimen_kk : 0.5*dimen_kk ) ; dset = EDIT_empty_copy( NULL ) ; EDIT_dset_items( dset , ADN_nxyz , iv_nxyz , ADN_xyzdel , fv_xyzdel , ADN_xyzorg , fv_xyzorg , ADN_xyzorient , iv_xyzorient , ADN_prefix , prefix , ADN_datum_all , datum , ADN_nvals , 1 , ADN_ntt , 0 , ADN_type , HEAD_FUNC_TYPE , ADN_func_type , FUNC_FIM_TYPE , ADN_none ) ; } if( THD_deathcon() && THD_is_file(DSET_HEADNAME(dset)) ) ERROR_exit("Output dataset already exists -- can't overwrite") ; if (do_head_only) { DSET_write_header(dset); exit(0); } /*-- make empty brick array for dataset --*/ EDIT_substitute_brick( dset , 0 , datum , NULL ) ; /* will make array */ nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxyz = nx*ny*nz; /* 19 Feb 2004: check and make mask if desired */ if( maskset != NULL && ( DSET_NX(maskset) != nx || DSET_NY(maskset) != ny || DSET_NZ(maskset) != nz ) ) ERROR_exit("-mask dataset doesn't match dimension of output dataset") ; if( maskset != NULL ){ mmask = THD_makemask( maskset , 0 , 1.0,-1.0 ) ; if( mmask == NULL ){ WARNING_message("Can't create mask for some reason!") ; } else { int nmask = THD_countmask( nxyz , mmask ) ; if( nmask == 0 ){ WARNING_message("0 voxels in mask -- ignoring it!") ; free((void *)mmask) ; mmask = NULL ; } else { INFO_message("%d voxels found in mask",nmask) ; } } DSET_delete(maskset) ; } /*-- fill new dataset brick with the -fval value --*/ switch( datum ){ case MRI_short: sbr = (short *) DSET_BRICK_ARRAY(dset,0) ; for( ii=0 ; ii < nxyz ; ii++ ) sbr[ii] = fval_short ; break ; case MRI_float: fbr = (float *) DSET_BRICK_ARRAY(dset,0) ; for( ii=0 ; ii < nxyz ; ii++ ) fbr[ii] = fval_float ; break ; case MRI_byte: bbr = (byte *) DSET_BRICK_ARRAY(dset,0) ; for( ii=0 ; ii < nxyz ; ii++ ) bbr[ii] = fval_byte ; break ; } /* 24 Nov 2000: get the bounding box for the dataset */ dx = fabs(dset->daxes->xxdel) ; if( dx <= 0.0 ) dx = 1.0 ; dy = fabs(dset->daxes->yydel) ; if( dy <= 0.0 ) dy = 1.0 ; dz = fabs(dset->daxes->zzdel) ; if( dz <= 0.0 ) dz = 1.0 ; #ifndef EXTEND_BBOX xxdown = dset->daxes->xxmin - 0.501 * dx ; xxup = dset->daxes->xxmax + 0.501 * dx ; yydown = dset->daxes->yymin - 0.501 * dy ; yyup = dset->daxes->yymax + 0.501 * dy ; zzdown = dset->daxes->zzmin - 0.501 * dz ; zzup = dset->daxes->zzmax + 0.501 * dz ; #else xxdown = dset->daxes->xxmin ; xxup = dset->daxes->xxmax ; yydown = dset->daxes->yymin ; yyup = dset->daxes->yymax ; zzdown = dset->daxes->zzmin ; zzup = dset->daxes->zzmax ; #endif /*-- loop over input files and read them line by line. --*/ /*(( iarg is already correct at the start of this loop ))*/ if( iarg == argc ) WARNING_message("No input files ==> creating empty dataset") ; if( ROImask != NULL ) INFO_message("Starting to fill via -ROImask dataset") ; else if( do_ijk ) INFO_message("Starting to fill via -ijk indexes") ; else INFO_message("Starting to fill via -xyz coordinates") ; for( ROInum=0 ; iarg < argc ; iarg++ ){ /*--- loop over files ---*/ /* get input file ready to read */ if( strcmp(argv[iarg],"-") == 0 || strncmp(argv[iarg],"stdin",5) == 0 ){ /* stdin */ fp = stdin ; } else { /* OK, open the damn file */ fp = fopen( argv[iarg] , "r" ) ; if( fp == NULL ){ WARNING_message("Can't open input file %s -- skipping it" , argv[iarg]) ; continue ; /* skip to end of iarg loop */ } } /* read lines, process and store */ ll = 0 ; while(1){ /*--- loop over lines ---*/ ll++ ; /* line count */ cp = afni_fgets( linbuf , NBUF , fp ) ; 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 ) continue ; /* all blanks */ if( linbuf[ii] == '/' ) continue ; /* comment */ if( linbuf[ii] == '#' ) continue ; /* comment */ if( linbuf[ii] == '<' ) continue ; /* XML? */ if( isalpha(linbuf[ii]) ) continue ; /* alphabetic */ /* changes commas to blanks [10 Nov 2008] */ for( jj=ii ; jj < kk ; jj++ ) if( linbuf[jj] == ',' || linbuf[jj] == ';' || linbuf[jj] == ':' ) linbuf[jj] = ' ' ; /* scan line for data */ vv = dval_float ; /* if not scanned in below, use the default value */ sv = dval_short ; bv = dval_byte ; vrad = srad ; /* 19 Feb 2004: default sphere radius */ nn = sscanf(linbuf+ii , "%f%f%f%f%f" , &xx,&yy,&zz,&vv,&vrad ) ; if( nn < 3 && ROImask == NULL ){ WARNING_message("File %s line %d: incomplete [%d]-- skipping",argv[iarg],ll,nn) ; continue ; } if( nn > 3 ){ if( MAYBE_SCAN(1,&vv) ){ WARNING_message("File %s line %d: replaced illegal value with 0", argv[iarg],ll ) ; } else if( isnan(vv) ){ WARNING_message("File %s line %d: keeping value NaN",argv[iarg],ll) ; } else { sv = SHORTIZE(vv) ; bv = BYTEIZE(vv) ; } if( thd_floatscan(1,&vrad) ){ WARNING_message("File %s line %d: replaced illegal radius with 0", argv[iarg],ll ) ; } } if( ROImask != NULL ){ /* 09 Nov 2011 -- new way of storing value */ int qq , nq ; ROInum++ ; for( nq=qq=0 ; qq < nxyz ; qq++ ){ /* loop over voxels */ if( ROIsss[qq] == ROInum ){ nq++ ; switch( datum ){ case MRI_float: fbr[qq] = xx ; vv = xx ; break ; case MRI_short: sv = SHORTIZE(xx) ; sbr[qq] = sv ; vv = sv ; break ; case MRI_byte: bv = BYTEIZE(xx) ; bbr[qq] = bv ; vv = bv ; break ; } } } ININFO_message(" filled %d voxel%s from ROI=%d with value=%g", nq , (nq==1) ? "\0" : "s" , ROInum , vv ) ; nfill += nq ; continue ; /* skip to next input line */ } /* get voxel index into (ii,jj,kk) */ if( do_ijk ){ /* inputs are (ii,jj,kk) themselves */ ii = (int) rint(xx) ; jj = (int) rint(yy) ; kk = (int) rint(zz) ; if( ii < 0 || ii >= nx ){ WARNING_message("File %s line %d: i index=%d is invalid", argv[iarg],ll,ii) ; continue ; } if( jj < 0 || jj >= ny ){ WARNING_message("File %s line %d: j index=%d is invalid", argv[iarg],ll,jj) ; continue ; } if( kk < 0 || kk >= nz ){ WARNING_message("File %s line %d: k index=%d is invalid", argv[iarg],ll,kk) ; continue ; } } else { /* inputs are coordinates => must convert to index */ THD_fvec3 mv , dv ; /* temp vectors */ THD_ivec3 iv ; THD_coorder_to_dicom( &cord , &xx,&yy,&zz ) ; /* to Dicom order */ LOAD_FVEC3( dv , xx,yy,zz ) ; mv = THD_dicomm_to_3dmm( dset , dv ) ; /* to Dataset order */ /* 24 Nov 2000: check (xx,yy,zz) for being inside the box */ if( mv.xyz[0] < xxdown || mv.xyz[0] > xxup ){ WARNING_message("File %s line %d: x coord=%g is outside %g .. %g" , argv[iarg],ll,mv.xyz[0] , xxdown,xxup ) ; continue ; } if( mv.xyz[1] < yydown || mv.xyz[1] > yyup ){ WARNING_message("File %s line %d: y coord=%g is outside %g .. %g" , argv[iarg],ll,mv.xyz[1] , yydown , yyup ) ; continue ; } if( mv.xyz[2] < zzdown || mv.xyz[2] > zzup ){ WARNING_message("File %s line %d: z coord=%g is outside %g .. %g" , argv[iarg],ll,mv.xyz[2] , zzdown , zzup ) ; continue ; } iv = THD_3dmm_to_3dind( dset , mv ) ; /* to Dataset index */ ii = iv.ijk[0]; jj = iv.ijk[1]; kk = iv.ijk[2]; /* save */ } /* now load individual voxel (ii,jj,kk) */ ijk = ii + jj*nx + kk*nx*ny ; if( mmask == NULL || mmask[ijk] ){ switch( datum ){ case MRI_float:{ if( fbr[ijk] == fval_float) nfill++ ; if( fbr[ijk] != fval_float && fbr[ijk] != vv ) WARNING_message("Overwrite float voxel %d %d %d",ii,jj,kk) ; fbr[ijk] = vv ; } break ; case MRI_short:{ if( sbr[ijk] == fval_short ) nfill++ ; if( sbr[ijk] != fval_short && sbr[ijk] != sv ) WARNING_message("Overwrite short voxel %d %d %d",ii,jj,kk) ; sbr[ijk] = sv ; } break ; case MRI_byte:{ if( bbr[ijk] == fval_byte ) nfill++ ; if( bbr[ijk] != fval_byte && bbr[ijk] != bv ) WARNING_message("Overwrite byte voxel %d %d %d",ii,jj,kk) ; bbr[ijk] = bv ; } break ; } } /* 19 Feb 2004: Make up radius of ellipsoid in voxel indexes Will put value vv into all voxels (aa,bb,cc) with (aa-ii)**2/qii + (bb-jj)**2/qjj + (cc-kk)**2/qkk <= 1.0 */ vrad *= 1.00001 ; /* expand a little */ rii = vrad/dx ; rjj = vrad/dy ; rkk = vrad/dz ; qii = rii*rii ; qjj = rjj*rjj ; qkk = rkk*rkk ; if( rii >= 1.0 || rjj >= 1.0 || rkk >= 1.0 ){ int aa,bb,cc , abot,atop,bbot,btop,cbot,ctop; float rr; byte assign_vox; assign_vox = 1; /* assume assigning for cubes */ abot = ii-(int)rint(rii) ; atop = ii+(int)rint(rii) ; if( abot < 0 ) abot = 0 ; if( atop >= nx ) atop = nx-1 ; bbot = jj-(int)rint(rjj) ; btop = jj+(int)rint(rjj) ; if( bbot < 0 ) bbot = 0 ; if( btop >= ny ) btop = ny-1 ; cbot = kk-(int)rint(rkk) ; ctop = kk+(int)rint(rkk) ; if( cbot < 0 ) cbot = 0 ; if( ctop >= nz ) ctop = nz-1 ; for( cc=cbot ; cc <= ctop ; cc++ ){ for( bb=bbot ; bb <= btop ; bb++ ){ for( aa=abot ; aa <= atop ; aa++ ){ if(!cubes) { /* cubes or spheres*/ rr = (aa-ii)*(aa-ii)/qii + (bb-jj)*(bb-jj)/qjj + (cc-kk)*(cc-kk)/qkk ; if( rr <= 1.00001 ) assign_vox = 1; else assign_vox = 0; } if( assign_vox ){ ijk = aa + bb*nx + cc*nx*ny ; /* (aa,bb,cc) in dataset */ if( mmask == NULL || mmask[ijk] ){ switch( datum ){ case MRI_float: if( fbr[ijk] == fval_float ) nfill++ ; fbr[ijk] = vv ; break ; case MRI_short: if( sbr[ijk] == fval_short ) nfill++ ; sbr[ijk] = sv ; break ; case MRI_byte: if( bbr[ijk] == fval_byte ) nfill++ ; bbr[ijk] = bv ; break ; } } } } } } } /* 19 Feb 2004: end of inserting a sphere */ } /*--- end of loop over input lines ---*/ /* close input file */ if( fp != stdin ) fclose( fp ) ; } /*--- end of loop over input files ---*/ if( !allow_NaN ) dset_floatscan(dset) ; tross_Make_History( "3dUndump" , argc,argv , dset ) ; INFO_message("Total number of voxels filled = %d",nfill) ; DSET_write(dset) ; INFO_message("Wrote out dataset %s",DSET_BRIKNAME(dset)) ; exit(0) ; }