#include "mrilib.h" /*------- Adapted from 3dTstat.c --------*/ static char prefix[THD_MAX_PREFIX] = "tsort" ; static int datum = MRI_float ; static int inc = 1 ; static void SORTS_tsfunc( double tzero , double tdelta , int npts , float ts[] , double ts_mean , double ts_slope , void *ud , int nbriks, float *val ) ; static void SORTS_itsfunc( double tzero , double tdelta , int npts , float ts[] , double ts_mean , double ts_slope , void *ud , int nbriks, float *val ) ; extern int *z_iqsort (float *x , int nx ); int main( int argc , char *argv[] ) { THD_3dim_dataset *old_dset , *new_dset ; /* input and output datasets */ int nopt, ii , nvals , rank; /*----- Help the pitiful user? -----*/ if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ printf("Usage: 3dTsort [options] dataset\n" "Sorts each voxel and produces a new dataset.\n" "\n" "Options:\n" " -prefix p = use string 'p' for the prefix of the\n" " output dataset [DEFAULT = 'tsort']\n" " -inc = sort into increasing order [default]\n" " -dec = sort into decreasing order\n" " -rank = output rank instead of sorted values\n" " ranks range from 1 to Nvals\n" " -val = output sorted values (default)\n" " -datum D = Coerce the output data to be stored as \n" " the given type D, which may be \n" " byte, short, or float (default). \n" "\n" "Notes:\n" "* Each voxel is sorted separately.\n" "* Sub-brick labels are not rearranged.\n" "* This program is useful only in limited cases.\n" " It was written to sort the -stim_times_IM\n" " beta weights output by 3dDeconvolve.\n" "* Also see program 1dTsort.\n" "\n" "Examples:\n" "setenv AFNI_1D_TIME YES\n" "echo '8 6 3 9 2 7' > test.1D\n" " 3dTsort -overwrite test.1D \n" " 1dcat tsort.1D\n" "\n" " 3dTsort -overwrite -rank test.1D \n" " 1dcat tsort.1D\n" "\n" " 3dTsort -overwrite -dec test.1D \n" " 1dcat tsort.1D\n" "\n" ) ; PRINT_COMPILE_DATE ; exit(0) ; } /* bureaucracy */ mainENTRY("3dTsort main"); machdep(); AFNI_logger("3dTsort",argc,argv); PRINT_VERSION("3dTsort"); AUTHOR("RW Cox"); /*--- scan command line for options ---*/ nopt = 1 ; rank = 0; while( nopt < argc && argv[nopt][0] == '-' ){ /*-- prefix --*/ if( strcmp(argv[nopt],"-prefix") == 0 ){ if( ++nopt >= argc ) ERROR_exit("-prefix needs an argument!\n"); MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ; if( !THD_filename_ok(prefix) ) ERROR_exit("%s is not a valid prefix!\n",prefix); nopt++ ; continue ; } /*-- -inc or -dec --*/ if( strncmp(argv[nopt],"-inc",4) == 0 ){ inc = 1 ; nopt++ ; continue ; } if( strncmp(argv[nopt],"-dec",4) == 0 ){ inc = 0 ; nopt++ ; continue ; } if( strncmp(argv[nopt],"-rank",5) == 0){ rank = 1; nopt++; continue; } if( strncmp(argv[nopt],"-val",4) == 0){ rank = 0; nopt++; continue; } if( strncasecmp(argv[nopt],"-datum",6) == 0 ){ if( ++nopt >= argc ) ERROR_exit("need an argument after -datum!\n") ; if( strcasecmp(argv[nopt],"short") == 0 ){ datum = MRI_short ; } else if( strcasecmp(argv[nopt],"float") == 0 ){ datum = MRI_float ; } else if( strcasecmp(argv[nopt],"byte") == 0 ){ datum = MRI_byte ; } else { ERROR_exit( "-datum of type '%s' not supported " "in 3dTsort!\n",argv[nopt]) ; } nopt++; continue; } /*-- Quien sabe'? --*/ ERROR_exit("Unknown option: %s\n",argv[nopt]) ; } /*----- read input dataset -----*/ if( nopt >= argc ) ERROR_exit(" No input dataset!?") ; old_dset = THD_open_dataset( argv[nopt] ) ; if( !ISVALID_DSET(old_dset) ) ERROR_exit("Can't open dataset %s\n",argv[nopt]); DSET_load(old_dset) ; if( !DSET_LOADED(old_dset) ) ERROR_exit("Can't load dataset %s\n",argv[nopt]) ; if (DSET_NUM_TIMES(old_dset)<2) { ERROR_exit( "Need at least 2 time points in series.\n" "Have only %d\n" "If using mutli-column 1D files, use\n" " setenv AFNI_1D_TIME YES\n" , DSET_NUM_TIMES(old_dset)); } nopt++ ; if( nopt < argc ) WARNING_message("Trailing inputs on command line ignored: %s ...",argv[nopt]) ; nvals = DSET_NVALS(old_dset) ; if( nvals < 2 ) ERROR_exit("Can't use dataset with < 2 values per voxel!\n") ; /*------------- ready to compute new dataset -----------*/ if (!rank) { new_dset = MAKER_4D_to_typed_fbuc( old_dset , /* input dataset */ prefix , /* output prefix */ datum , /* output datum */ 0 , /* ignore count */ 0 , /* don't detrend */ nvals , /* number of briks */ SORTS_tsfunc , /* timeseries processor */ NULL /* data for tsfunc */ ) ; } else { new_dset = MAKER_4D_to_typed_fbuc( old_dset , /* input dataset */ prefix , /* output prefix */ datum , /* output datum */ 0 , /* ignore count */ 0 , /* don't detrend */ nvals , /* number of briks */ SORTS_itsfunc , /* timeseries processor */ NULL /* data for tsfunc */ ) ; } if( new_dset != NULL ){ tross_Copy_History( old_dset , new_dset ) ; tross_Make_History( "3dTsort" , argc,argv , new_dset ) ; if( DSET_NUM_TIMES(old_dset) > 1 ) EDIT_dset_items( new_dset , ADN_ntt , DSET_NVALS(old_dset) , ADN_ttorg , DSET_TIMEORIGIN(old_dset) , ADN_ttdel , DSET_TR(old_dset) , ADN_tunits , UNITS_SEC_TYPE , NULL ) ; DSET_write( new_dset ) ; WROTE_DSET( new_dset ) ; } else { ERROR_exit("Unable to compute output dataset!\n") ; } exit(0) ; } /********************************************************************** Function that does the real work ***********************************************************************/ static void SORTS_tsfunc( double tzero, double tdelta , int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val ) { int ii , nval ; /** is this a "notification"? **/ if( val == NULL ){ #if 0 if( npts > 0 ){ /* the "start notification" */ } else { /* the "end notification" */ } #endif return ; } /** do the work **/ nval = MIN(nbriks,npts) ; memcpy( val , ts , sizeof(float)*nval ) ; if( inc == 0 ){ for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ; } qsort_float( nval , val ) ; if( inc == 0 ){ for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ; } return ; } static void SORTS_itsfunc( double tzero, double tdelta , int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val ) { int ii , nval ; int *rnk; /** is this a "notification"? **/ if( val == NULL ){ #if 0 if( npts > 0 ){ /* the "start notification" */ } else { /* the "end notification" */ } #endif return ; } /** do the work **/ nval = MIN(nbriks,npts) ; memcpy( val , ts , sizeof(float)*nval ) ; /* Using an inverse sorting function, for excitement */ #if 0 /* the dumber way */ for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ; rnk = z_iqsort( val, nval ) ; if( inc == 1 ){ for( ii=0 ; ii < nval ; ii++ ) val[ii] = rnk[ii] +1 ; }else { for( ii=0 ; ii < nval ; ii++ ) val[ii] = nval - rnk[ii] ; } #else /* the dumb way */ rnk = z_iqsort( val, nval ) ; if( inc == 1 ){ for( ii=0 ; ii < nval ; ii++ ) val[rnk[ii]] = nval - ii ; }else { for( ii=0 ; ii < nval ; ii++ ) val[rnk[ii]] = ii+1 ; } #endif free(rnk); rnk=NULL; return ; }