#include "mrilib.h" /* [PT: Oct 15, 2018] + added '-extent_ijk_to_file FF' option to get slice numbers of auto-bboxing in a nice text file. [PT: Oct 15, 2018] + change where+how input dset check occurs, so subbrick selection is possible [PT: Oct 15, 2018] + added a couple other slice info things: '-extent_ijk' to put bounding box slices to screen '-extent_ijk_midslice' to put midslice of bounding box to screen [PT: Jan 17, 2023] + added a couple other slice info things: '-extent_xyz_to_file' to put bounding box XYZ coords to file + removed this note, because it does not appear to be the case (anymore?) " Also note that this value is calculated before\n" " any '-npad ...' option, so it would ignore that.\n" + make help file formatting a bit more readable [PT: Apr 18, 2024] + new opt for simpler scripting: -extent_xyz_quiet */ /* A function that takes a dataset as input, and based on its orientation outputs a string of 'i', 'j' and 'k' values. For example, the orientation 'RPI' -> 'ijk', and 'AIL' -> 'jki'. */ int THD_fill_char_int_3_ijk(THD_dataxes * daxes, char ostr[4]); int THD_fill_char_int_3_ijk(THD_dataxes * daxes, char ostr[4]) { int nn; char *ref_str = "ijk\0"; if ( ! daxes || ! ostr ) return 1; /* INFO_message("TEST1: %d %d %d", ORIENT_xyzint[daxes->xxorient]-1, ORIENT_xyzint[daxes->yyorient]-1, ORIENT_xyzint[daxes->zzorient]-1); INFO_message("TEST1: %c %c %c", ref_str[ORIENT_xyzint[daxes->xxorient]-1 ], ref_str[ORIENT_xyzint[daxes->yyorient]-1 ], ref_str[ORIENT_xyzint[daxes->zzorient]-1 ]); */ ostr[0] = ref_str[ORIENT_xyzint[daxes->xxorient]-1 ]; ostr[1] = ref_str[ORIENT_xyzint[daxes->yyorient]-1 ]; ostr[2] = ref_str[ORIENT_xyzint[daxes->zzorient]-1 ]; return 0; }; void help_autobox() { printf( "Usage: 3dAutobox [options] DATASET\n" "Computes size of a box that fits around the volume.\n" "Also can be used to crop the volume to that box.\n" "\n" "The default 'info message'-based terminal text is a set of IJK coords.\n" "See below for options to display coordinates in other ways, as well as\n" "to save them in a text file. Please note in particular the difference\n" "between *ijk* and *ijkord* outputs, for scripting.\n" "\n" "OPTIONS: ~1~\n" "\n" "-prefix PREFIX :Crop the input dataset to the size of the box, and\n" " write an output dataset with PREFIX for the name.\n" " *If -prefix is not used, no new volume is written out,\n" " just the (x,y,z) extents of the voxels to be kept.\n" "\n" "-input DATASET :An alternate way to specify the input dataset.\n" " The default method is to pass DATASET as\n" " the last parameter on the command line.\n" "\n" "-noclust :Don't do any clustering to find box. Any non-zero\n" " voxel will be preserved in the cropped volume.\n" " The default method uses some clustering to find the\n" " cropping box, and will clip off small isolated blobs.\n" "\n" "-extent :Write to standard out the spatial extent of the box\n" "\n" "-extent_xyz_quiet :The same numbers as '-extent', but only numbers and\n" " no string content. Ordering is RLAPIS.\n" "\n" "-extent_ijk :Write out the 6 auto bbox ijk slice numbers to\n" " screen:\n" " imin imax jmin jmax kmin kmax\n" " Note that resampling would affect the ijk vals (but\n" " not necessarily the xyz ones).\n" "\n" "-extent_ijk_to_file FF\n" " :Write out the 6 auto bbox ijk slice numbers to\n" " a simple-formatted text file FF (single row file):\n" " imin imax jmin jmax kmin kmax\n" " (same notes as above apply).\n" "\n" "-extent_ijk_midslice :Write out the 3 ijk midslices of the autobox to\n" " the screen:\n" " imid jmid kmid\n" " These are obtained via: (imin + imax)/2, etc.\n" "\n" "-extent_ijkord :Write out the 6 auto bbox ijk slice numbers to screen\n" " but in a particular order and format (see 'NOTE on\n" " *ijkord* format', below).\n" " NB: This ordering is useful if you want to use\n" " the output indices in 3dcalc expressions.\n" "\n" "-extent_ijkord_to_file FFORRD\n" " :Write out the 6 auto bbox ijk slice numbers to a file\n" " but in a particular order and format (see 'NOTE on\n" " *ijkord* format', below).\n" " NB: This option is quite useful if you want to use\n" " the output indices in 3dcalc expressions.\n" "\n" "-extent_xyz_to_file GG\n" " :Write out the 6 auto bbox xyz coords to\n" " a simple-formatted text file GG (single row file):\n" " xmin xmax ymin ymax zmin zmax\n" " (same values as '-extent').\n" "\n" "-extent_xyz_midslice :Write out the 3 xyz midslices of the autobox to\n" " the screen:\n" " xmid ymid zmid\n" " These are obtained via: (xmin + xmax)/2, etc.\n" " These follow the same meaning as '-extent'.\n" "\n" "-npad NNN :Number of extra voxels to pad on each side of box,\n" " since some troublesome people (that's you, LRF) want\n" " this feature for no apparent reason.\n" " ** With this option, it is possible to get a dataset\n" " thatis actually bigger than the input.\n" " ** You can input a negative value for NNN, which will\n" " crop the dataset even more than the automatic method.\n" "\n" "-npad_safety_on :Constrain npad-ded extents to be within dset. So, \n" " each index is bounded to be in range [0, L-1], where L\n" " is matrix length along that dimension.\n" "\n" "\n" "NOTE on *ijkord* format ~1~\n" "\n" "Using any of the '-*ijkord*' options above will output pairs of ijk\n" "indices just like the regular ijk options, **but** they will be ordered\n" "in a way that you can associate each of the i, j, and k indices with\n" "a standard x, y and z coordinate direction. Without this ordering,\n" "resampling a dataset could change what index is associated with which\n" "coordinate axis. That situation can be confusing for scripting (and\n" "by confusing, we mean 'bad').\n" "\n" "The output format for any '-*ijkord*' options is a 3x3 table, where\n" "the first column is the index value (i, j or k), and the next two\n" "columns are the min and max interval boundaries for the autobox.\n" "Importantly, the rows are placed in order so that the top corresponds\n" "to the x-axis, the middle to the y-axis and the bottom to the z-axis.\n" "\n" "So, if you had the following table output for a dset:\n" " k 10 170\n" " i 35 254\n" " j 21 199\n" "\n" "... you would look at the third row for the min/max slice values\n" "along the z-axis, and you would use the index 'j' to refer to it in,\n" "say, a 3dcalc expression.\n" "\n" "Note that the above example table output came from a dataset with ASL\n" "orientation. We can see how that fits, recalling that the first,\n" "second and third rows tell us about x, y and z info, respectively; and\n" "that i, j and k refer to the first, second and third characters in the\n" "orientation string. So, the third (z-like) row contains a j, which\n" "points us at the middle character in the orientation, which is S, which\n" "is along the z-axis---all consistent! Similarly, the top (x-like) row\n" "contains a k, which points us at the last char in the orientation,\n" "which is L and that is along the x-axis---phew!\n" "\n" "The main point of this would be to extra this information and use it\n" "in a script. If you knew that you wanted the z-slice range to use\n" "in a 3dcalc 'within()' expression, then you could extract the 3rd row\n" "to get the correct index and slice ranges, e.g., in tcsh:\n" " set vvv = `sed -n 3p FILE_ijkord.txt`\n" "\n" "... where now ${vvv} will have 3 values, the first of which is the\n" "relevant index letter, then the min and max slice range values.\n" "So an example 3dcalc expression to keep values only within\n" "that slice range:\n" " 3dcalc \\\n" " -a DSET \\\n" " -expr \"a*within(${vvv[1]},${vvv[2]},${vvv[3]})\" \\\n" " -prefix DSET_SUBSET\n" "\n" ) ; PRINT_COMPILE_DATE ; return ; } /*----------------------------------------------------------------------------*/ int main( int argc , char * argv[] ) { THD_3dim_dataset *dset, *outset=NULL; int iarg=1, npad = 0, extent=0; char *prefix=NULL, *iname=NULL; char *oijkext = NULL; FILE *fout_ijkext = NULL; char *oxyzext = NULL; FILE *fout_xyzext = NULL; char *oijkordext = NULL; FILE *fout_ijkordext = NULL; int extent_xyz_q=0; int extent_ijk=0; int extent_ijk_midslice=0; int extent_ijkord=0; int imid=0, jmid=0, kmid=0; int extent_xyz_midslice=0; float xmid=0, ymid=0, zmid=0; char *full_orient = NULL; /* PT: default output ordering of IJK depends on the dset orientation, so we don't actually which indices go along the z-direction, like if we wanted to make a 3dcalc expression using 'k' to select rows. This array contains the IJK extents ordered so the first pair is along the x-direction, second pair along the y-direction and third along z-direction. Used with the "-*ijkord*" opts */ int ijk_as_rai[6] = {-1, -1, -1, -1, -1}; char ijk_order[4] = "---\0"; /* PT: default npad can give indices outside dset. Activating this flag with the '-npad_safety_on' opt will constrain each index to be in range [0, L-1], where L is matrix size along that dim. */ int DO_NPAD_SAFELY = 0; /*-- startup bureaucracy --*/ mainENTRY("3dAutobox main"); machdep(); AFNI_logger("3dAutobox",argc,argv); PRINT_VERSION("3dAutobox") ; /*-- read command line options --*/ iarg = 1 ; while( iarg < argc && argv[iarg][0] == '-' ){ if( strcmp(argv[iarg],"-help") == 0 || strcmp(argv[iarg],"-h") == 0){ help_autobox(); exit(0) ; } if( strcmp(argv[iarg],"-prefix") == 0 ){ prefix = argv[++iarg] ; if( !THD_filename_ok(prefix) ){ fprintf(stderr,"** 3dAutobox: Illegal string after -prefix!\n"); exit(1) ; } iarg++ ; continue ; } if( strcmp(argv[iarg],"-input") == 0 ){ iname = argv[++iarg] ; // This is a bad check, because it doesn't permit subbrick // selection! Will do check later. /*if( !THD_filename_ok(iname) ){ fprintf(stderr,"** 3dAutobox: Illegal string after -input!\n"); exit(1) ; }*/ iarg++ ; continue ; } if( strcmp(argv[iarg],"-noclust") == 0 ){ MRI_autobbox_clust(0) ; /* turn of clustering and clipping */ THD_autobbox_clip(0) ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-npad") == 0 ){ npad = (int)strtod(argv[++iarg],NULL) ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-npad_safety_on") == 0 ){ DO_NPAD_SAFELY = 1; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_ijk") == 0 ){ extent_ijk = 1 ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_ijk_to_file") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after '-extent_ijk_to_file'\n") ; oijkext = argv[iarg]; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_ijkord") == 0 ){ extent_ijkord = 1 ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_ijkord_to_file") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after '-extent_ijkord_to_file'\n") ; oijkordext = argv[iarg]; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_xyz_to_file") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after '-extent_xyz_to_file'\n") ; oxyzext = argv[iarg]; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_xyz_quiet") == 0 ){ extent_xyz_q = 1; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_ijk_midslice") == 0 ){ extent_ijk_midslice = 1 ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent_xyz_midslice") == 0 ){ extent_xyz_midslice = 1 ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-extent") == 0 ){ extent = 1 ; iarg++ ; continue ; } /*- washappenin, dood? -*/ ERROR_message("** 3dAutobox: %s makes no sense here.\n", argv[iarg]) ; suggest_best_prog_option(argv[0], argv[iarg]); exit(1); } if( argc < 2){ help_autobox(); exit(0); } /* got input ? */ if( iarg == argc-1 ) iname = argv[iarg] ; else if (iarg != argc) ERROR_exit("** 3dAutobox: %s is nonsense on the line \n" " I know you're John; stop pretending you have an accent!", argv[iarg]) ; if( !iname ) ERROR_exit("** 3dAutobox: Where is my input?") ; /*-- read data --*/ dset = THD_open_dataset(iname); // Check here instead of after -input. if( dset == NULL ) ERROR_exit("Can't open time series dataset '%s'.",iname); CHECK_OPEN_ERROR(dset,iname); if( DSET_BRICK_TYPE(dset,0) != MRI_short && DSET_BRICK_TYPE(dset,0) != MRI_byte && DSET_BRICK_TYPE(dset,0) != MRI_float ) ERROR_exit("** ILLEGAL dataset type: %s :-(", MRI_type_name[DSET_BRICK_TYPE(dset,0)]) ; DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ; { int nx=DSET_NX(dset), ny=DSET_NY(dset), nz=DSET_NZ(dset), nxy=nx*ny ; int xm=-1,xp=-1,ym=-1,yp=-1,zm=-1,zp=-1; THD_autobbox( dset , &xm,&xp , &ym,&yp , &zm,&zp, NULL ) ; xm -= npad; ym -= npad; zm -= npad; /* for LRF */ xp += npad; yp += npad; zp += npad; /* Constrain each index to be in [0, L-1], where L is matr len. Some checks here might be extraneous but... they are all valid */ if ( DO_NPAD_SAFELY ) { if ( xm < 0 ) xm = 0; if ( xm >= nx ) xm = nx - 1; if ( xp < 0 ) xp = 0; if ( xp >= nx ) xp = nx - 1; if ( ym < 0 ) ym = 0; if ( ym >= ny ) ym = ny - 1; if ( yp < 0 ) yp = 0; if ( yp >= ny ) yp = ny - 1; if ( zm < 0 ) zm = 0; if ( zm >= nz ) zm = nz - 1; if ( zp < 0 ) zp = 0; if ( zp >= nz ) zp = nz - 1; } INFO_message("Auto bbox: x=%d..%d y=%d..%d z=%d..%d\n", xm,xp,ym,yp,zm,zp ) ; // [PT: Jan 17, 2023] do the work for ijkord, including getting // the associated string values { int i; int axes_ord[3] = {-10, -10, -10}; char test_str[] = "---\0"; // get int orient codes i = THD_fill_orient_int_3_rlpais(dset->daxes, axes_ord); // ... and convert the codes, ordering for x (=0), y (=1) and z (=2) for( i=0 ; i<3 ; i++ ) axes_ord[i] = ORIENT_xyzint[axes_ord[i]] - 1; // place the paired ijk values into xyz-like order ijk_as_rai[2*axes_ord[0]] = xm; ijk_as_rai[2*axes_ord[0]+1] = xp; ijk_as_rai[2*axes_ord[1]] = ym; ijk_as_rai[2*axes_ord[1]+1] = yp; ijk_as_rai[2*axes_ord[2]] = zm; ijk_as_rai[2*axes_ord[2]+1] = zp; THD_fill_char_int_3_ijk(dset->daxes, test_str); //INFO_message("TEST: %s", test_str); for( i=0 ; i<3 ; i++ ) ijk_order[i] = test_str[axes_ord[i]]; //INFO_message("TESTb: %s", ijk_order); //INFO_message("Auto bbox: iord=%d..%d jord=%d..%d kord=%d..%d\n", // ijk_as_rai[0], ijk_as_rai[1], ijk_as_rai[2], // ijk_as_rai[3], ijk_as_rai[4], ijk_as_rai[5] ); } // [PT: Oct 18, 2018] New output text file, if desired if( oijkext ) { if( (fout_ijkext = fopen(oijkext, "w")) == NULL ) { fprintf(stderr, "\n\nError opening file %s.\n", oijkext); exit(1); } fprintf( fout_ijkext, "%8d %8d %8d %8d %8d %8d\n", xm, xp, ym, yp, zm, zp ); fclose(fout_ijkext); INFO_message("Wrote ijk extents file: %s", oijkext); } if( extent_ijk ) printf( "%8d %8d %8d %8d %8d %8d\n", xm, xp, ym, yp, zm, zp ); if( extent_ijk_midslice ) { imid = (xm + xp) / 2; // integer division fine, b/c we need ints jmid = (ym + yp) / 2; kmid = (zm + zp) / 2; printf( "%8d %8d %8d\n", imid, jmid, kmid ); } if( extent_ijkord ) printf( "%c %8d %8d\n%c %8d %8d\n%c %8d %8d\n", ijk_order[0], ijk_as_rai[0], ijk_as_rai[1], ijk_order[1], ijk_as_rai[2], ijk_as_rai[3], ijk_order[2], ijk_as_rai[4], ijk_as_rai[5] ); if ( (extent && !prefix) || (extent_xyz_midslice && !prefix) || \ (oxyzext && !prefix) || (oijkordext && !prefix) || \ (extent_xyz_q && !prefix) ) prefix = "EXTENT_ONLY"; if( prefix ){ outset = THD_zeropad( dset , -xm, xp-nx+1, -ym, yp-ny+1, -zm, zp-nz+1, prefix , ZPAD_IJK ) ; if( THD_deathcon() && THD_is_file(DSET_HEADNAME(outset)) ) ERROR_exit("3dAutobox: output file %s already exists :-(", DSET_HEADNAME(outset) ) ; if( outset == NULL ) ERROR_exit("3dAutobox: Some error occurred in processing :-(") ; tross_Copy_History( dset , outset ) ; /* 31 Jan 2001 - RWCox */ tross_Make_History( "3dAutobox" , argc,argv , outset ) ; if (!strstr(prefix,"EXTENT_ONLY")) { DSET_write(outset) ; INFO_message("3dAutobox: output dataset = %s", DSET_BRIKNAME(outset)) ; } if (extent) { float RL_AP_IS[6]; THD_dset_extent(outset, '-', RL_AP_IS); printf("Extent auto bbox: R=%f L=%f A=%f P=%f I=%f S=%f\n", RL_AP_IS[0],RL_AP_IS[1], RL_AP_IS[2],RL_AP_IS[3], RL_AP_IS[4],RL_AP_IS[5] ) ; } // [PT: April 18, 2024] only 6 numbers output RLAPIS if( extent_xyz_q ) { float RL_AP_IS[6]; THD_dset_extent(outset, '-', RL_AP_IS); printf( "%f %f %f %f %f %f\n", RL_AP_IS[0],RL_AP_IS[1], RL_AP_IS[2],RL_AP_IS[3], RL_AP_IS[4],RL_AP_IS[5] ) ; } // [PT: Jan 17, 2023] New output text file: xyz if( oxyzext ) { float RL_AP_IS[6]; THD_dset_extent(outset, '-', RL_AP_IS); if( (fout_xyzext = fopen(oxyzext, "w")) == NULL ) { fprintf(stderr, "\n\nError opening file %s.\n", oxyzext); exit(1); } fprintf( fout_xyzext, "%f %f %f %f %f %f\n", RL_AP_IS[0], RL_AP_IS[1], RL_AP_IS[2], RL_AP_IS[3], RL_AP_IS[4], RL_AP_IS[5] ); fclose(fout_xyzext); INFO_message("Wrote xyz extents file: %s", oxyzext); } // [PT: Jan 17, 2023] New output text file: ijkord if( oijkordext ) { if( (fout_ijkordext = fopen(oijkordext, "w")) == NULL ) { fprintf(stderr, "\n\nError opening file %s.\n", oijkordext); exit(1); } fprintf( fout_ijkordext, "%c %8d %8d\n%c %8d %8d\n%c %8d %8d\n", ijk_order[0], ijk_as_rai[0], ijk_as_rai[1], ijk_order[1], ijk_as_rai[2], ijk_as_rai[3], ijk_order[2], ijk_as_rai[4], ijk_as_rai[5] ); fclose(fout_ijkordext); INFO_message("Wrote ijkord extents file: %s", oijkordext); } if( extent_xyz_midslice ) { float RL_AP_IS2[6]; THD_dset_extent(outset, '-', RL_AP_IS2); xmid = (RL_AP_IS2[0] + RL_AP_IS2[1]) / 2.; ymid = (RL_AP_IS2[2] + RL_AP_IS2[3]) / 2.; zmid = (RL_AP_IS2[4] + RL_AP_IS2[5]) / 2.; printf( "%10.5f %10.5f %10.5f\n", xmid, ymid, zmid ); } } } exit(0) ; }