#include "mrilib.h" int main( int argc , char * argv[] ) { float mfrac=0.50f , val ; MRI_IMAGE *medim ; THD_3dim_dataset *dset ; char *gprefix=NULL ; int iarg=1 , doall=0 ; if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ printf( "Usage: 3dClipLevel [options] dataset\n" "Estimates the value at which to clip the anatomical dataset so\n" " that background regions are set to zero.\n" "\n" "The program's output is a single number sent to stdout. This\n" " value can be 'captured' to a shell variable using the backward\n" " single quote operator; a trivial csh/tcsh example is\n" "\n" " set ccc = `3dClipLevel -mfrac 0.333 Elvis+orig`\n" " 3dcalc -a Elvis+orig -expr \"step(a-$ccc)\" -prefix Presley\n" "\n" "Algorithm:\n" " (a) Set some initial clip value using wizardry (AKA 'variance').\n" " (b) Find the median of all positive values >= clip value.\n" " (c) Set the clip value to 0.50 of this median.\n" " (d) Loop back to (b) until the clip value doesn't change.\n" "This method was made up out of nothing, based on histogram gazing.\n" "\n" "Options:\n" "--------\n" " -mfrac ff = Use the number ff instead of 0.50 in the algorithm.\n" " -doall = Apply the algorithm to each sub-brick separately.\n" " [Cannot be combined with '-grad'!]\n" "\n" " -grad ppp = In addition to using the 'one size fits all routine',\n" " also compute a 'gradual' clip level as a function\n" " of voxel position, and output that to a dataset with\n" " prefix 'ppp'.\n" " [This is the same 'gradual' clip level that is now the\n" " default in 3dAutomask - as of 24 Oct 2006.\n" " You can use this option to see how 3dAutomask clips\n" " the dataset as its first step. The algorithm above is\n" " is used in each octant of the dataset, and then these\n" " 8 values are interpolated to cover the whole volume.]\n" "Notes:\n" "------\n" "* Use at your own risk! You might want to use the AFNI Histogram\n" " plugin to see if the results are reasonable. This program is\n" " likely to produce bad results on images gathered with local\n" " RF coils, or with pulse sequences with unusual contrasts.\n" "\n" "* For brain images, most brain voxels seem to be in the range from\n" " the clip level (mfrac=0.5) to about 3-3.5 times the clip level.\n" " - In T1-weighted images, voxels above that level are usually\n" " blood vessels (e.g., inflow artifact brightens them).\n" "\n" "* If the input dataset has more than 1 sub-brick, the data is\n" " analyzed on the median volume -- at each voxel, the median\n" " of all sub-bricks at that voxel is computed, and then this\n" " median volume is used in the histogram algorithm.\n" "\n" "* If the input dataset is short- or byte-valued, the output will\n" " be an integer; otherwise, the output is a float value.\n" "\n" "* Example -- Scaling a sequence of sub-bricks from a collection of\n" " anatomicals from different sites to have about the\n" " same numerical range (from 0 to 255):\n" " 3dTcat -prefix input anat_*+tlrc.HEAD\n" " 3dClipLevel -doall input+tlrc > clip.1D\n" " 3dcalc -datum byte -nscale -a input+tlrc -b clip.1D \\\n" " -expr '255*max(0,min(1,a/(3.2*b)))' -verb -prefix scaled\n" "----------------------------------------------------------------------\n" "* Author: Emperor Zhark -- Sadistic Galactic Domination since 1994!\n\n" ) ; PRINT_COMPILE_DATE ; exit(0) ; } mainENTRY("3dClipLevel main"); machdep(); AFNI_logger("3dCliplevel",argc,argv); /*-- options --*/ while( iarg < argc && argv[iarg][0] == '-' ){ if( strcmp(argv[iarg],"-doall") == 0 ){ /* 01 Nov 2010 */ doall++ ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-mfrac") == 0 || strcmp(argv[iarg],"-clfrac") == 0 ){ mfrac = strtod( argv[++iarg] , NULL ) ; if( mfrac <= 0.0f ) ERROR_exit("Illegal -mfrac '%s'",argv[iarg]) ; while( mfrac >= 1.0f ) mfrac *= 0.01f ; iarg++ ; continue ; } if( strcmp(argv[iarg],"-grad") == 0 ){ gprefix = argv[++iarg] ; if( !THD_filename_ok(gprefix) ) ERROR_exit("Illegal -grad '%s'",argv[iarg]) ; iarg++ ; continue ; } ERROR_exit("Unknown option: %s\n",argv[iarg]) ; } /*-- read input dataset --*/ dset = THD_open_dataset(argv[iarg]) ; CHECK_OPEN_ERROR(dset,argv[iarg]) ; DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ; if( !doall ){ /*-- get median at each voxel --*/ medim = THD_median_brick( dset ) ; if( medim == NULL ) ERROR_exit("Can't load dataset '%s'",argv[iarg]); DSET_unload(dset) ; /* no longer needed */ val = THD_cliplevel( medim , mfrac ) ; /* floating point clip level */ if( !THD_need_brick_factor(dset) && /* convert to integer? */ DSET_datum_constant(dset) && (DSET_BRICK_TYPE(dset,0)==MRI_short || DSET_BRICK_TYPE(dset,0)==MRI_byte) ) val = (float)rint((double)val) ; printf("%g\n",val) ; /***** write the output value! *****/ /*--- 25 Oct 2006: create the gradual clip dataset? ---*/ if( gprefix != NULL && val > 0.0f ){ THD_3dim_dataset *gset ; MRI_IMAGE *gim ; gim = THD_cliplevel_gradual( medim , mfrac ) ; if( gim == NULL ) ERROR_exit("Can't compute gradual clip?!") ; gset = EDIT_empty_copy(dset) ; EDIT_dset_items( gset , ADN_nvals , 1 , ADN_ntt , 0 , ADN_prefix, gprefix , ADN_none ) ; EDIT_substitute_brick( gset , 0 , MRI_float , MRI_FLOAT_PTR(gim) ) ; tross_Copy_History( dset , gset ) ; tross_Make_History( "3dClipLevel" , argc,argv , gset ) ; DSET_write(gset) ; DSET_delete(gset) ; } } else { /* -doall [01 Nov 2010] */ int nvals=DSET_NVALS(dset) , ii ; MRI_IMAGE *bim ; printf("# %s %d values\n",DSET_BRIKNAME(dset),nvals) ; for( ii=0 ; ii < nvals ; ii++ ){ bim = THD_extract_float_brick(ii,dset) ; DSET_unload_one(dset,ii) ; val = THD_cliplevel( bim , mfrac ) ; /* floating point clip level */ printf("%g\n",val) ; mri_free(bim) ; } } exit(0) ; }