1 /*********************** 3dRank.c **********************************************/
2 /* Author: Ziad Saad May 2009 */
3 #include "mrilib.h"
4 #include "uthash.h"
5
6 typedef struct {
7 int id; /* keep it named 'id' to facilitate use of convenience
8 macros in uthash . */
9 UT_hash_handle hh; /* keep name for same reason */
10 int index;
11 } INT_HASH_DATUM;
12
13
14 extern int * UniqueInt (int *y, int ysz, int *kunq, int Sorted );
15
Rank_help(void)16 void Rank_help(void) {
17 printf(
18 "Usage: 3dRank [-prefix PREFIX] <-input DATASET1 [DATASET2 ...]>\n"
19 "\n"
20 "Replaces voxel values by their rank in the set of\n"
21 "values collected over all voxels in all input datasets\n"
22 "If you input one dataset, the output should be identical\n"
23 "to the -1rank option in 3dmerge\n"
24 "\n"
25 "This program only works on datasets of integral storage type, \n"
26 "and on integral valued data stored as floats.\n"
27 "\n"
28 " -input DATASET1 [DATASET2 ...]: Input datasets.\n"
29 " Acceptable data types are:\n"
30 " byte, short, and floats.\n"
31 " -prefix PREFIX: Ouput prefix.\n"
32 " If you have multiple datasets on input\n"
33 " the prefix is preceded by r00., r01., etc.\n"
34 " If no prefix is given, the default is \n"
35 " rank.DATASET1, rank.DATASET2, etc.\n"
36 "\n"
37 " In addition to the ranked volume, a rank map\n"
38 " 1D file is created. It shows the mapping from \n"
39 " the rank (1st column) to the integral values \n"
40 " (2nd column) in the input dataset. Sub-brick float \n"
41 " factors are ignored.\n"
42 "\n"
43 " -ver = print author and version info\n"
44 " -help = print this help screen\n"
45 ) ;
46 printf("\n" MASTER_SHORTHELP_STRING ) ;
47 PRINT_COMPILE_DATE ;
48 return;
49 }
50
51
52
53 /*! Replace a voxel's value by the value's rank in the entire set of input datasets */
main(int argc,char * argv[])54 int main( int argc , char * argv[] )
55 {
56 THD_3dim_dataset ** dsets_in = NULL, *dset=NULL; /*input and output datasets*/
57 int nopt=0, nbriks=0, nsubbriks=0, ib=0, isb=0;
58 byte *cmask=NULL;
59 int *all_uniques=NULL, **uniques=NULL, *final_unq=NULL, *N_uniques=NULL;
60 int N_final_unq=0, iun=0, total_unq=0;
61 INT_HASH_DATUM *rmap=NULL, *hd=NULL;
62 int imax=0, iunq=0, ii=0, id = 0;
63 long int off=0;
64 char *prefix=NULL;
65 char stmp[THD_MAX_PREFIX+1]={""};
66 FILE *fout=NULL;
67
68 /*----- Read command line -----*/
69 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
70 Rank_help ();
71 exit(0) ;
72 }
73
74 mainENTRY("3dRank main"); machdep(); AFNI_logger("3dRank",argc,argv);
75 nopt = 1 ;
76
77 while( nopt < argc && argv[nopt][0] == '-' ){
78 if( strcmp(argv[nopt],"-ver") == 0 ){
79 PRINT_VERSION("3dRank"); AUTHOR("Ziad Saad");
80 nopt++; continue;
81 }
82
83 if( strcmp(argv[nopt],"-help") == 0 ){
84 Rank_help();
85 exit(0) ;
86 }
87
88 if( strcmp(argv[nopt],"-prefix") == 0 ){
89 ++nopt;
90 if (nopt>=argc) {
91 fprintf(stderr,"**ERROR: Need string after -prefix\n");
92 exit(1);
93 }
94 prefix = argv[nopt] ;
95 ++nopt; continue;
96 }
97 if( strcmp(argv[nopt],"-input") == 0 ){
98 dsets_in = (THD_3dim_dataset**)
99 calloc(argc-nopt+1, sizeof(THD_3dim_dataset*));
100 ++nopt; nbriks=0;
101 while (nopt < argc ) {
102 dsets_in[nbriks] = THD_open_dataset( argv[nopt] );
103 if( !ISVALID_DSET(dsets_in[nbriks]) ){
104 fprintf(stderr,"**ERROR: can't open dataset %s\n",argv[nopt]) ;
105 exit(1);
106 }
107 ++nopt; ++nbriks;
108 }
109 continue;
110 }
111
112 ERROR_exit( " Error - unknown option %s", argv[nopt]);
113 }
114 if (nopt < argc) {
115 ERROR_exit( " Error unexplained trailing option: %s\n", argv[nopt]);
116 }
117 if (!nbriks) {
118 ERROR_exit( " Error no volumes entered on command line?");
119 }
120
121 /* some checks and inits*/
122 nsubbriks = 0;
123 for (ib = 0; ib<nbriks; ++ib) {
124 if (!is_integral_dset(dsets_in[ib], 0)) {
125 ERROR_exit("Dset %s is not of an integral data type.",
126 DSET_PREFIX(dsets_in[ib]));
127 }
128 nsubbriks += DSET_NVALS(dsets_in[ib]);
129 }
130
131 /* Now get unique arrays */
132 uniques = (int **)calloc(nsubbriks, sizeof(int*));
133 N_uniques = (int *)calloc(nsubbriks, sizeof(int));
134 total_unq = 0;
135 iun = 0;
136 for (ib = 0; ib<nbriks; ++ib) {
137 DSET_mallocize(dsets_in[ib]); DSET_load(dsets_in[ib]);
138 for (isb=0; isb<DSET_NVALS(dsets_in[ib]); ++isb) {
139 uniques[iun] = THD_unique_vals(dsets_in[ib], isb,
140 &(N_uniques[iun]), cmask);
141 total_unq += N_uniques[iun];
142 ++iun;
143 }
144 }
145
146 /* put all the arrays together and get the unique of the uniques */
147 all_uniques = (int *)calloc(total_unq, sizeof(int));
148 off=0;
149 for (iun=0; iun<nsubbriks; ++iun) {
150 memcpy(all_uniques+off, uniques[iun], N_uniques[iun]*sizeof(int));
151 off += N_uniques[iun];
152 }
153
154 /* free intermediate unique arrays */
155 for (iun=0; iun<nsubbriks; ++iun) {
156 free(uniques[iun]);
157 }
158 free(uniques); uniques=NULL;
159 free(N_uniques); N_uniques=NULL;
160
161 /* get unique of catenated array */
162 if (!(final_unq = UniqueInt (all_uniques, total_unq, &N_final_unq, 0 ))) {
163 ERROR_exit( " Failed to get unique list (%d, %d, %d) ",
164 total_unq, N_final_unq, nsubbriks);
165 }
166 free(all_uniques); all_uniques=NULL;
167
168 if (prefix) {
169 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
170 "%s.rankmap.1D", prefix);
171 } else {
172 if (nbriks == 1) {
173 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
174 "%s.rankmap.1D", DSET_PREFIX(dsets_in[0]));
175 } else {
176 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
177 "%s+.rankmap.1D", DSET_PREFIX(dsets_in[0]));
178 }
179 }
180
181 if (stmp[0]) {
182 if ((fout = fopen(stmp,"w"))) {
183 fprintf(fout, "#Rank Map (%d unique values)\n", N_final_unq);
184 fprintf(fout, "#Col. 0: Rank\n");
185 fprintf(fout, "#Col. 1: Input Dset Value\n");
186 }
187 }
188
189
190 /* get the maximum integer in the unique array */
191 imax = 0;
192 for (iunq=0; iunq<N_final_unq; ++iunq) {
193 if (final_unq[iunq] > imax) imax = final_unq[iunq];
194 if (fout) fprintf(fout, "%d %d\n", iunq, final_unq[iunq]);
195 hd = (INT_HASH_DATUM*)calloc(1,sizeof(INT_HASH_DATUM));
196 hd->id = final_unq[iunq];
197 hd->index = iunq;
198 HASH_ADD_INT(rmap, id, hd);
199 }
200
201 fclose(fout); fout=NULL;
202
203 /* now cycle over all dsets and replace their voxel values with rank */
204 for (ib = 0; ib<nbriks; ++ib) {
205 for (isb=0; isb<DSET_NVALS(dsets_in[ib]); ++isb) {
206 EDIT_BRICK_LABEL( dsets_in[ib],isb, "rank" ) ;
207 EDIT_BRICK_TO_NOSTAT( dsets_in[ib],isb ) ;
208 EDIT_BRICK_FACTOR( dsets_in[ib],isb, 0.0);/* no factors for rank*/
209 switch (DSET_BRICK_TYPE(dsets_in[ib],isb) ){
210 default:
211 fprintf(stderr,
212 "** Bad dset type for unique operation.\n"
213 "Only Byte, Short, and float dsets are allowed.\n");
214 break ; /* this should not happen here,
215 so don't bother returning*/
216 case MRI_short:{
217 short *mar = (short *) DSET_ARRAY(dsets_in[ib],isb) ;
218 if (imax > MRI_TYPE_maxval[MRI_short]) {
219 WARNING_message("Maximum rank value of %d is more\n"
220 "than maximum value for dset datatype of %d\n",
221 imax, MRI_TYPE_maxval[MRI_short]);
222 }
223 for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )
224 if (!cmask || cmask[ii]) {
225 id = (int)mar[ii];
226 HASH_FIND_INT(rmap,&id ,hd);
227 if (hd) mar[ii] = (short)(hd->index);
228 else
229 ERROR_exit("** Failed to find key %d in hash table\n",id);
230 } else mar[ii] = 0;
231 }
232 break ;
233 case MRI_byte:{
234 byte *mar ;
235 if (imax > MRI_TYPE_maxval[MRI_byte]) {
236 WARNING_message("Maximum rank value of %d is more\n"
237 "than maximum value for dset datatype of %d\n",
238 imax, MRI_TYPE_maxval[MRI_byte]);
239 }
240 mar = (byte *) DSET_ARRAY(dsets_in[ib],isb) ;
241 for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )
242 if (!cmask || cmask[ii]) {
243 id = (int)mar[ii];
244 HASH_FIND_INT(rmap,&id ,hd);
245 if (hd) mar[ii] = (byte)(hd->index);
246 else
247 ERROR_exit("** Failed to find key %d in hash table\n",id);
248 } else mar[ii] = 0;
249 }
250 break ;
251 case MRI_float:{
252 float *mar = (float *) DSET_ARRAY(dsets_in[ib],isb) ;
253 for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )
254 if (!cmask || cmask[ii]) {
255 id = (int)mar[ii]; /* Assuming float is integral valued */
256 HASH_FIND_INT(rmap,&id ,hd);
257 if (hd) mar[ii] = (float)(hd->index);
258 else
259 ERROR_exit("** Failed to find key %d in hash table\n",id);
260 } else mar[ii] = 0;
261 }
262 break ;
263
264 }
265 }
266
267 /* update range, etc. */
268 THD_load_statistics(dsets_in[ib]);
269
270 /* Now write the bricks */
271 if (prefix) {
272 if (nbriks == 1) {
273 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
274 "%s", prefix);
275 } else {
276 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
277 "r%02d.%s", ib, prefix);
278 }
279 } else {
280 snprintf(stmp, sizeof(char)*THD_MAX_PREFIX,
281 "rank.%s", DSET_PREFIX(dsets_in[ib]));
282 }
283
284 EDIT_dset_items( dsets_in[ib] ,
285 ADN_prefix , stmp ,
286 ADN_none ) ;
287
288 /* storage_mode was set in EDIT_dset_items(), based on the prefix,
289 do not reset it here [21 Nov 2019 rickr] */
290 /* dsets_in[ib]->dblk->diskptr->storage_mode = STORAGE_BY_BRICK; */
291
292 tross_Make_History( "3dRank" , argc, argv , dsets_in[ib] ) ;
293 if (DSET_IS_MASTERED(dsets_in[ib])) {
294 /* THD_write_3dim_dataset won't write a mastered dude */
295 dset = EDIT_full_copy(dsets_in[ib],stmp);
296 } else {
297 dset = dsets_in[ib];
298 }
299
300 /* New ID */
301 ZERO_IDCODE(dset->idcode);
302 dset->idcode = MCW_new_idcode() ;
303
304 if (!THD_write_3dim_dataset( NULL, stmp, dset,True )) {
305 ERROR_message("Failed to write %s", stmp);
306 exit(1);
307 } else {
308 WROTE_DSET(dsets_in[ib]);
309 if (dset != dsets_in[ib]) DSET_deletepp(dset);
310 DSET_deletepp(dsets_in[ib]);
311
312 }
313 }
314
315 /* destroy hash */
316 while (rmap) {
317 hd = rmap;
318 HASH_DEL(rmap,hd);
319 free(hd);
320 }
321
322 free(final_unq); final_unq=NULL;
323
324 exit(0);
325 }
326
327