1 #include "mrilib.h"
Help_3dRankizer(TFORM targ,int detail)2 int Help_3dRankizer(TFORM targ, int detail)
3 {
4
5 if (detail >= 0) {
6 sphinx_printf(targ,
7 "Usage: 3dRankizer [options] dataset\n"
8 "Output = Rank of each voxel as sorted into increasing value.\n"
9 " - Ties get the average rank.\n"
10 " - Not the same as 3dRank!\n"
11 " - Only sub-brick #0 is processed at this time!\n"
12 " - Ranks start at 1 and increase:\n"
13 " Input = 0 3 4 4 7 9\n"
14 " Output = 1 2 3.5 3.5 5 6\n"
15 "Options:\n"
16 " -brank bbb Set the 'base' rank to 'bbb' instead of 1.\n"
17 " (You could also do this with 3dcalc.)\n"
18 " -mask mset Means to use the dataset 'mset' as a mask:\n"
19 " Only voxels with nonzero values in 'mset'\n"
20 " will be used from 'dataset'. Voxels outside\n"
21 " the mask will get rank 0.\n"
22 " -prefix ppp Write results into float-format dataset 'ppp'\n"
23 " Output is in float format to allow for\n"
24 " non-integer ranks resulting from ties.\n"
25 " -percentize : Divide rank by the number of voxels in the dataset x 100.0 \n"
26 " -percentize_mask : Divide rank by the number of voxels in the mask x 100.0 \n"
27 "\n"
28 "Author: RW Cox [[a quick hack for his own purposes]]\n"
29 );
30 }
31 if (detail >= 1) {
32 PRINT_COMPILE_DATE ;
33 }
34 return(0);
35 }
36
main(int argc,char * argv[])37 int main( int argc , char * argv[] )
38 {
39 THD_3dim_dataset *inset=NULL , *outset=NULL , *mask_dset=NULL ;
40 MRI_IMAGE *fim=NULL ; float *far=NULL ;
41 int iarg , ndset , nvox , ii , mcount=0 , percentize = 0;
42 char *prefix = "rankizer" ;
43 byte *mmm=NULL ;
44 float brank=1.0f ;
45
46
47 /*---- official startup ---*/
48
49 PRINT_VERSION("3dRankizer"); mainENTRY("3dRankizer main"); machdep();
50 AFNI_logger("3dRankizer",argc,argv); AUTHOR("Zhark of the Ineffable Rank");
51
52 /*-- read command line arguments --*/
53 /*-- command line scan --*/
54
55 iarg = 1 ;
56 while( iarg < argc && argv[iarg][0] == '-' ){
57 CHECK_HELP(arg,Help_3dRankizer);
58
59 if( strncmp(argv[iarg],"-brank",5) == 0 ){
60 if( iarg+1 >= argc )
61 ERROR_exit("-brank option requires a following argument!") ;
62 brank = (float)strtod(argv[++iarg],NULL) ;
63 iarg++ ; continue ;
64 }
65
66 if( strncmp(argv[iarg],"-mask",5) == 0 ){
67 if( mask_dset != NULL )
68 ERROR_exit("Cannot have two -mask options!") ;
69 if( iarg+1 >= argc )
70 ERROR_exit("-mask option requires a following argument!") ;
71 mask_dset = THD_open_dataset( argv[++iarg] ) ;
72 if( mask_dset == NULL )
73 ERROR_exit("Cannot open mask dataset!") ;
74 iarg++ ; continue ;
75 }
76
77 if( strcmp(argv[iarg],"-prefix") == 0 ){
78 if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix'") ;
79 prefix = strdup(argv[iarg]) ;
80 if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal value after '-prefix'") ;
81 iarg++ ; continue ;
82 }
83
84 if( strcmp(argv[iarg],"-percentize") == 0 ) {
85 percentize = 1;
86 iarg++ ; continue ;
87 }
88
89 if( strcmp(argv[iarg],"-percentize_mask") == 0 ) {
90 percentize = 2;
91 iarg++ ; continue ;
92 }
93
94 ERROR_message("Unknown option: %s\n",argv[iarg]) ;
95 suggest_best_prog_option(argv[0], argv[iarg]);
96 exit(1);
97 }
98
99 if( argc < 3 ){
100 Help_3dRankizer(TXT, 0);
101 PRINT_COMPILE_DATE ; exit(0) ;
102 }
103
104 /* should have 1 more arg */
105
106 ndset = argc - iarg ;
107 if( ndset < 1 ) ERROR_exit("No input dataset!?") ;
108 else if( ndset > 1 ) WARNING_message("Too many input datasets!") ;
109
110 inset = THD_open_dataset( argv[iarg] ) ; CHECK_OPEN_ERROR(inset,argv[iarg]) ;
111 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
112 fim = THD_extract_float_brick(0,inset) ; DSET_unload(inset) ;
113 far = MRI_FLOAT_PTR(fim) ;
114 nvox= DSET_NVOX(inset) ;
115
116 /* make a byte mask from mask dataset */
117
118 if( mask_dset != NULL ){
119 if( DSET_NVOX(mask_dset) != nvox )
120 ERROR_exit("Input and mask datasets are not same dimensions!");
121 mmm = THD_makemask( mask_dset , 0 , 1.0f,-1.0f ) ;
122 mcount = THD_countmask( nvox , mmm ) ;
123 INFO_message("%d voxels in the mask",mcount) ;
124 if( mcount <= 5 ) ERROR_exit("Mask is too small!") ;
125 DSET_delete(mask_dset) ;
126 }
127
128 if( mmm == NULL ){
129 mcount = nvox;
130 rank_order_float( nvox , far ) ;
131 for( ii=0 ; ii < nvox ; ii++ ) far[ii] += brank ;
132 } else {
133 float fmin=far[0] ;
134 for( ii=1 ; ii < nvox ; ii++ ) if( far[ii] < fmin ) fmin = far[ii] ;
135 if( fmin > 0.0f ) fmin = 0.0f ;
136 else if( fmin == 0.0f ) fmin = -1.0f ;
137 else fmin = -2.0f*fmin-1.0f ;
138 for( ii=0 ; ii < nvox ; ii++ ) if( !mmm[ii] ) far[ii] = fmin ;
139 rank_order_float( nvox , far ) ;
140 fmin = (nvox-mcount) - brank ;
141 for( ii=0 ; ii < nvox ; ii++ ){
142 if( mmm[ii] ) far[ii] = far[ii] - fmin ;
143 else far[ii] = 0.0f ;
144 }
145 }
146
147 if (percentize) {
148 float fac=1.0;
149 if (percentize == 1) fac = 100.0/nvox;
150 else if (percentize == 2) fac = 100.0/mcount;
151 for( ii=0 ; ii < nvox ; ii++ ) far[ii] *= fac;
152 }
153
154 outset = EDIT_empty_copy( inset ) ;
155 EDIT_dset_items( outset ,
156 ADN_prefix , prefix ,
157 ADN_brick_fac , NULL ,
158 ADN_nvals , 1 ,
159 ADN_ntt , 0 ,
160 ADN_none ) ;
161 EDIT_substitute_brick( outset , 0 , MRI_float , far ) ;
162 DSET_write(outset) ; WROTE_DSET(outset) ;
163 exit(0) ;
164 }
165