1 #include "mrilib.h"
2
3 /*------- Adapted from 3dTstat.c --------*/
4
5 static char prefix[THD_MAX_PREFIX] = "tsort" ;
6 static int datum = MRI_float ;
7 static int inc = 1 ;
8 static int do_random = 0 ;
9
10 static void SORTS_tsfunc( double tzero , double tdelta ,
11 int npts , float ts[] , double ts_mean ,
12 double ts_slope , void *ud , int nbriks, float *val ) ;
13 static void SORTS_itsfunc( double tzero , double tdelta ,
14 int npts , float ts[] , double ts_mean ,
15 double ts_slope , void *ud , int nbriks, float *val ) ;
16 static void SORTS_rtsfunc( double tzero , double tdelta ,
17 int npts , float ts[] , double ts_mean ,
18 double ts_slope , void *ud , int nbriks, float *val ) ;
19 static void SORTS_ranfunc( double tzero , double tdelta ,
20 int npts , float ts[] , double ts_mean ,
21 double ts_slope , void *ud , int nbriks, float *val ) ;
22 static void SORTS_FFTfunc( double tzero , double tdelta ,
23 int npts , float ts[] , double ts_mean ,
24 double ts_slope , void *ud , int nbriks, float *val ) ;
25 extern int *z_iqsort (float *x , int nx );
26
27 static unsigned short xran_pm[3] = { 23456 , 34567 , 54321 } ;
28 #define SET_XRAN(xr,rss) \
29 ( (xr)[0]=((rss)>>16), (xr)[1]=((rss)&65535), (xr)[2]=(xr)[0]+(xr)[1]+17 )
30
main(int argc,char * argv[])31 int main( int argc , char *argv[] )
32 {
33 THD_3dim_dataset *old_dset , *new_dset ; /* input and output datasets */
34 int nopt, ii , nvals , rank=0 ;
35
36 /*----- Help the pitiful user? -----*/
37
38 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
39 printf("Usage: 3dTsort [options] dataset\n"
40 "Sorts each voxel and produces a new dataset.\n"
41 "\n"
42 "Options:\n"
43 " -prefix p = use string 'p' for the prefix of the\n"
44 " output dataset [DEFAULT = 'tsort']\n"
45 " -inc = sort into increasing order [default]\n"
46 " -dec = sort into decreasing order\n"
47 " -rank = output rank instead of sorted values\n"
48 " ranks range from 1 to Nvals\n"
49 " -ind = output sorting index. (0 to Nvals -1)\n"
50 " See example below.\n"
51 " -val = output sorted values (default)\n"
52 " -random = randomly shuffle (permute) the time points in each voxel\n"
53 " * Each voxel is permuted independently!\n"
54 " * Why is this here? Someone asked for it :)\n"
55 " -ranFFT = randomize each time series by scrambling the FFT phase\n"
56 " * Each voxel is treated separately!\n"
57 " * Why is this here? cf. Matthew 7:7-8 :)\n"
58 " -ranDFT = Almost the same as above, but:\n"
59 " * In '-ranFFT', the FFT length is taken\n"
60 " to be the next integer >= data length\n"
61 " for which the FFT algorithm is efficient.\n"
62 " This will result in data padding unless\n"
63 " the data length is exactly 'nice' for FFT.\n"
64 " * In '-ranDFT', the DFT length is exactly\n"
65 " the data length. If the data length is\n"
66 " a large-ish prime number (say 997), this\n"
67 " operation can be slow.\n"
68 " * The DFT/FFT algorithm is reasonably fast\n"
69 " when the data length prime factors contain\n"
70 " only 2s, 3s, and/or 5s.\n"
71 " * Using '-ranDFT' can preserve the spectral\n"
72 " (temporal correlation) structure of the\n"
73 " original data a little better than '-ranFFT'.\n"
74 " * The only reason to use '-ranFFT' instead of\n"
75 " '-ranDFT' is for speed. For example, with\n"
76 " 997 time points, '-ranFFT' was about 13 times\n"
77 " faster (FFT length=1000) than '-ranDFT'.\n"
78 " -datum D = Coerce the output data to be stored as \n"
79 " the given type D, which may be \n"
80 " byte, short, or float (default). \n"
81 "\n"
82 "Notes:\n"
83 "* Each voxel is sorted (or processed) separately.\n"
84 "* Sub-brick labels are not rearranged!\n"
85 "* This program is useful only in limited cases.\n"
86 " It was written to sort the -stim_times_IM\n"
87 " beta weights output by 3dDeconvolve.\n"
88 "* Also see program 1dTsort, for sorting text files of numbers.\n"
89 "\n"
90 "Examples:\n"
91 "setenv AFNI_1D_TIME YES\n"
92 "echo '8 6 3 9 2 7' > test.1D\n"
93 " 3dTsort -overwrite test.1D \n"
94 " 1dcat tsort.1D\n"
95 "\n"
96 " 3dTsort -overwrite -rank test.1D \n"
97 " 1dcat tsort.1D\n"
98 "\n"
99 "\n"
100 " 3dTsort -overwrite -ind test.1D \n"
101 " 1dcat tsort.1D\n"
102 "\n"
103 " 3dTsort -overwrite -dec test.1D \n"
104 " 1dcat tsort.1D\n"
105 "\n"
106 ) ;
107 PRINT_COMPILE_DATE ; exit(0) ;
108 }
109
110 /* bureaucracy */
111
112 mainENTRY("3dTsort main"); machdep(); AFNI_logger("3dTsort",argc,argv);
113 PRINT_VERSION("3dTsort"); AUTHOR("RW Cox");
114
115 /*--- scan command line for options ---*/
116
117 nopt = 1 ;
118 rank = 0;
119 while( nopt < argc && argv[nopt][0] == '-' ){
120
121 /*-- prefix --*/
122
123 if( strcmp(argv[nopt],"-prefix") == 0 ){
124 if( ++nopt >= argc ) ERROR_exit("-prefix needs an argument!\n");
125 MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
126 if( !THD_filename_ok(prefix) )
127 ERROR_exit("%s is not a valid prefix!\n",prefix);
128 nopt++ ; continue ;
129 }
130
131 /*-- -inc or -dec --*/
132
133 if( strncmp(argv[nopt],"-inc",4) == 0 ){
134 inc = 1 ; nopt++ ; continue ;
135 }
136 if( strncmp(argv[nopt],"-dec",4) == 0 ){
137 inc = 0 ; nopt++ ; continue ;
138 }
139 if( strncmp(argv[nopt],"-rank",5) == 0){
140 rank = 1; nopt++; continue;
141 }
142 if( strncmp(argv[nopt],"-val",4) == 0){
143 rank = 0; nopt++; continue;
144 }
145 if( strncmp(argv[nopt],"-ind",5) == 0){
146 rank = -1; nopt++; continue;
147 }
148 if( strncmp(argv[nopt],"-random",6) == 0){
149 unsigned int spm ;
150 rank = 0; do_random = 1;
151 init_rand_seed(0) ;
152 spm = lrand48() ;
153 SET_XRAN(xran_pm,spm) ;
154 nopt++; continue;
155 }
156 if( strncmp(argv[nopt],"-ranFFT",6) == 0){
157 unsigned int spm ;
158 rank = 0; do_random = -1;
159 init_rand_seed(0) ;
160 spm = lrand48() ;
161 SET_XRAN(xran_pm,spm) ;
162 nopt++; continue;
163 }
164 if( strncmp(argv[nopt],"-ranDFT",6) == 0){
165 unsigned int spm ;
166 rank = 0; do_random = -2;
167 init_rand_seed(0) ;
168 spm = lrand48() ;
169 SET_XRAN(xran_pm,spm) ;
170 nopt++; continue;
171 }
172 if( strncasecmp(argv[nopt],"-datum",6) == 0 ){
173 if( ++nopt >= argc )
174 ERROR_exit("need an argument after -datum!\n") ;
175 if( strcasecmp(argv[nopt],"short") == 0 ){
176 datum = MRI_short ;
177 } else if( strcasecmp(argv[nopt],"float") == 0 ){
178 datum = MRI_float ;
179 } else if( strcasecmp(argv[nopt],"byte") == 0 ){
180 datum = MRI_byte ;
181 } else {
182 ERROR_exit( "-datum of type '%s' not supported "
183 "in 3dTsort!\n",argv[nopt]) ;
184 }
185 nopt++; continue;
186 }
187 /*-- Quien sabe'? --*/
188
189 ERROR_exit("Unknown option: %s\n",argv[nopt]) ;
190 }
191
192 /*----- read input dataset -----*/
193
194 if( nopt >= argc ) ERROR_exit(" No input dataset!?") ;
195
196 old_dset = THD_open_dataset( argv[nopt] ) ;
197 if( !ISVALID_DSET(old_dset) )
198 ERROR_exit("Can't open dataset %s\n",argv[nopt]);
199 DSET_load(old_dset) ;
200 if( !DSET_LOADED(old_dset) )
201 ERROR_exit("Can't load dataset %s\n",argv[nopt]) ;
202 if (DSET_NUM_TIMES(old_dset)<2) {
203 ERROR_exit( "Need at least 2 time points in series.\n"
204 "Have only %d\n"
205 "If using mutli-column 1D files, use\n"
206 " setenv AFNI_1D_TIME YES\n"
207 , DSET_NUM_TIMES(old_dset));
208 }
209 nopt++ ;
210 if( nopt < argc )
211 WARNING_message("Trailing inputs on command line ignored: %s ...",argv[nopt]) ;
212
213 nvals = DSET_NVALS(old_dset) ;
214 if( nvals < 2 )
215 ERROR_exit("Can't use dataset with < 2 values per voxel!\n") ;
216
217 /*------------- ready to compute new dataset -----------*/
218
219 if( do_random == 1 ){ /* shuffling */
220 new_dset = MAKER_4D_to_typed_fbuc(
221 old_dset , /* input dataset */
222 prefix , /* output prefix */
223 datum , /* output datum */
224 0 , /* ignore count */
225 0 , /* don't detrend */
226 nvals , /* number of briks */
227 SORTS_ranfunc , /* timeseries processor */
228 NULL, /* data for tsfunc */
229 NULL, /* mask */
230 0 /* Allow auto scaling of output */
231 ) ;
232 } else if ( do_random < 0 ){ /* scrambling */
233 new_dset = MAKER_4D_to_typed_fbuc(
234 old_dset , /* input dataset */
235 prefix , /* output prefix */
236 datum , /* output datum */
237 0 , /* ignore count */
238 0 , /* don't detrend */
239 nvals , /* number of briks */
240 SORTS_FFTfunc , /* timeseries processor */
241 NULL, /* data for tsfunc */
242 NULL, /* mask */
243 0 /* Allow auto scaling of output */
244 ) ;
245 } else if (!rank) {
246 new_dset = MAKER_4D_to_typed_fbuc(
247 old_dset , /* input dataset */
248 prefix , /* output prefix */
249 datum , /* output datum */
250 0 , /* ignore count */
251 0 , /* don't detrend */
252 nvals , /* number of briks */
253 SORTS_tsfunc , /* timeseries processor */
254 NULL, /* data for tsfunc */
255 NULL, /* mask */
256 0 /* Allow auto scaling of output */
257 ) ;
258 } else if (rank == 1) {
259 new_dset = MAKER_4D_to_typed_fbuc(
260 old_dset , /* input dataset */
261 prefix , /* output prefix */
262 datum , /* output datum */
263 0 , /* ignore count */
264 0 , /* don't detrend */
265 nvals , /* number of briks */
266 SORTS_itsfunc , /* timeseries processor */
267 NULL, /* data for tsfunc */
268 NULL, /* mask */
269 0 /* Allow auto scaling of output */
270 ) ;
271 } else {
272 new_dset = MAKER_4D_to_typed_fbuc(
273 old_dset , /* input dataset */
274 prefix , /* output prefix */
275 datum , /* output datum */
276 0 , /* ignore count */
277 0 , /* don't detrend */
278 nvals , /* number of briks */
279 SORTS_rtsfunc , /* timeseries processor */
280 NULL, /* data for tsfunc */
281 NULL, /* mask */
282 0 /* Allow auto scaling of output */
283 ) ;
284 }
285 if( new_dset != NULL ){
286 tross_Copy_History( old_dset , new_dset ) ;
287 tross_Make_History( "3dTsort" , argc,argv , new_dset ) ;
288 if( DSET_NUM_TIMES(old_dset) > 1 )
289 EDIT_dset_items( new_dset ,
290 ADN_ntt , DSET_NVALS(old_dset) ,
291 ADN_ttorg , DSET_TIMEORIGIN(old_dset) ,
292 ADN_ttdel , DSET_TR(old_dset) ,
293 ADN_tunits , UNITS_SEC_TYPE ,
294 NULL ) ;
295 DSET_write( new_dset ) ; WROTE_DSET( new_dset ) ;
296 } else {
297 ERROR_exit("Unable to compute output dataset!\n") ;
298 }
299
300 exit(0) ;
301 }
302
303 /**********************************************************************
304 Function that does the real work
305 ***********************************************************************/
306
SORTS_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)307 static void SORTS_tsfunc( double tzero, double tdelta ,
308 int npts, float ts[],
309 double ts_mean, double ts_slope,
310 void *ud, int nbriks, float *val )
311 {
312 int ii , nval ;
313
314 /** is this a "notification"? **/
315
316 if( val == NULL ){
317 #if 0
318 if( npts > 0 ){ /* the "start notification" */
319 } else { /* the "end notification" */
320 }
321 #endif
322 return ;
323 }
324
325 /** do the work **/
326
327 nval = MIN(nbriks,npts) ;
328 memcpy( val , ts , sizeof(float)*nval ) ;
329 if( inc == 0 ){
330 for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ;
331 }
332 qsort_float( nval , val ) ;
333 if( inc == 0 ){
334 for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ;
335 }
336 return ;
337 }
338
339 /*--------------------------------------------------------------------*/
340
SORTS_itsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)341 static void SORTS_itsfunc( double tzero, double tdelta ,
342 int npts, float ts[],
343 double ts_mean, double ts_slope,
344 void *ud, int nbriks, float *val )
345 {
346 int ii , nval ;
347 int *rnk;
348 /** is this a "notification"? **/
349
350 if( val == NULL ){
351 #if 0
352 if( npts > 0 ){ /* the "start notification" */
353 } else { /* the "end notification" */
354 }
355 #endif
356 return ;
357 }
358
359 /** do the work **/
360
361 nval = MIN(nbriks,npts) ;
362 memcpy( val , ts , sizeof(float)*nval ) ;
363
364 /* Using an inverse sorting function, for excitement */
365 #if 0 /* the dumber way */
366 for( ii=0 ; ii < nval ; ii++ ) val[ii] = -val[ii] ;
367 rnk = z_iqsort( val, nval ) ;
368 if( inc == 1 ){
369 for( ii=0 ; ii < nval ; ii++ ) val[ii] = rnk[ii] +1 ;
370 }else {
371 for( ii=0 ; ii < nval ; ii++ ) val[ii] = nval - rnk[ii] ;
372 }
373 #else /* the dumb way */
374 rnk = z_iqsort( val, nval ) ;
375 if( inc == 1 ){
376 for( ii=0 ; ii < nval ; ii++ ) val[rnk[ii]] = nval - ii ;
377 }else {
378 for( ii=0 ; ii < nval ; ii++ ) val[rnk[ii]] = ii+1 ;
379 }
380 #endif
381 free(rnk); rnk=NULL;
382 return ;
383 }
384
385 /*--------------------------------------------------------------------*/
386
SORTS_rtsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)387 static void SORTS_rtsfunc( double tzero, double tdelta ,
388 int npts, float ts[],
389 double ts_mean, double ts_slope,
390 void *ud, int nbriks, float *val )
391 {
392 int ii , nval ;
393 int *rnk;
394 /** is this a "notification"? **/
395
396 if( val == NULL ){
397 return ;
398 }
399
400 /** do the work **/
401
402 nval = MIN(nbriks,npts) ;
403 memcpy( val , ts , sizeof(float)*nval ) ;
404
405 /* Using an inverse sorting function, for excitement */
406 rnk = z_iqsort( val, nval ) ;
407 if( inc == 1 ){
408 for( ii=0 ; ii < nval ; ii++ ) val[nval - ii-1] = rnk[ii];
409 }else {
410 for( ii=0 ; ii < nval ; ii++ ) val[ii] = rnk[ii] ;
411 }
412 free(rnk); rnk=NULL;
413 return ;
414 }
415
416 /*--------------------------------------------------------------------*/
417
418 static int p_nxy = 0 ; /* total length of data */
419 static float *p_xyar = NULL ; /* array to hold both samples */
420 static int *p_ijar = NULL ; /* permutation array */
421
float_permute(int nxy,float * ar)422 static void float_permute( int nxy , float *ar )
423 {
424 int ii,jj,tt ;
425
426 if( nxy < 2 ) return ; /* how did this happen? */
427
428 if( nxy > p_nxy ){ /* make workspaces */
429 p_nxy = nxy ;
430 p_xyar = (float *)realloc(p_xyar,sizeof(float)*p_nxy) ;
431 p_ijar = (int *)realloc(p_ijar,sizeof(int )*p_nxy) ;
432 }
433
434 /* initialize the permutation a little randomly */
435
436 tt = nrand48(xran_pm) % p_nxy ;
437 for( ii=0 ; ii < p_nxy ; ii++ ) p_ijar[ii] = (ii+tt)%p_nxy ;
438
439 /* create a random-ish permutation */
440 /* https://en.wikipedia.org/wiki/Random_permutation */
441
442 for( ii=0 ; ii < p_nxy-1 ; ii++ ){
443 jj = (nrand48(xran_pm)>>3) % (p_nxy-ii) ; /* jj in 0..p_nxy-ii-1 inclusive */
444 /* so ii+jj in ii..p_nxy-1 inclusive */
445 if( jj > 0 ){ /* swap */
446 tt = p_ijar[ii] ; p_ijar[ii] = p_ijar[ii+jj] ; p_ijar[ii+jj] = tt ;
447 }
448 }
449
450 for( ii=0 ; ii < p_nxy ; ii++ ) p_xyar[ii] = ar[ii] ;
451 for( ii=0 ; ii < p_nxy ; ii++ ) ar[ii] = p_xyar[p_ijar[ii]] ;
452
453 return ;
454 }
455
SORTS_ranfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)456 static void SORTS_ranfunc( double tzero, double tdelta ,
457 int npts, float ts[],
458 double ts_mean, double ts_slope,
459 void *ud, int nbriks, float *val )
460 {
461 int ii , nval ;
462
463 /** is this a "notification"? **/
464
465 if( val == NULL ){
466 #if 0
467 if( npts > 0 ){ /* the "start notification" */
468 } else { /* the "end notification" */
469 }
470 #endif
471 return ;
472 }
473
474 /** do the work **/
475
476 nval = MIN(nbriks,npts) ;
477 memcpy( val , ts , sizeof(float)*nval ) ;
478 float_permute( nval , val ) ;
479 return ;
480 }
481
482 /*--------------------------------------------------------------------*/
483
484 static int s_nx=0 , s_nfft=0 , s_nfft2=0 , s_nftop=0 ;
485 static complex *s_cxar=NULL ;
486 static float **s_pref=NULL ;
487 static float *s_ffar=NULL ;
488 static float *s_ffit=NULL ;
489
490 #define CXROT(cc,th) \
491 do{ complex rr; \
492 rr.r = cos(th); rr.i = sin(th); \
493 CMULTBY(cc,rr) ; } while(0)
494
495 #define TPI 6.283185f
496
float_scramble(int nx,float * ar)497 static void float_scramble( int nx , float *ar )
498 {
499 int ii,jj ; float theta ;
500
501 if( nx < 5 ) return ; /* how did this happen? */
502
503 if( nx != s_nx ){
504 s_nx = nx ;
505 if( do_random == -1 ){
506 s_nfft = csfft_nextup_even(s_nx) ; s_nfft2 = s_nfft/2 ;
507 } else {
508 s_nfft = s_nx ;
509 s_nfft2 = (s_nfft %2 == 0 ) ? s_nfft/2 : 0;
510 }
511 s_nftop = s_nfft / 2 ;
512 s_cxar = (complex *)realloc(s_cxar,sizeof(complex)*s_nfft) ;
513 s_ffar = (float * )realloc(s_ffar,sizeof(float) *s_nfft) ;
514 if( s_pref != NULL ){
515 free(s_pref[0]) ; free(s_pref[1]) ; free(s_pref[2]) ; free(s_pref) ;
516 }
517 s_pref = THD_build_polyref( 3 , s_nx ) ;
518 if( s_ffit == NULL ) s_ffit = (float *)malloc(sizeof(float)*3) ;
519 csfft_scale_inverse(1) ;
520 INFO_message("time series length = %d ==> FFT length = %d",
521 s_nx,s_nfft) ;
522 }
523
524 memcpy( s_ffar , ar , sizeof(float)*s_nx ) ;
525
526 THD_generic_detrend_LSQ( s_nx , s_ffar , -1 , 3 , s_pref , s_ffit ) ;
527
528 for( ii=0 ; ii < s_nfft ; ii++ ){
529 s_cxar[ii].r = s_ffar[ii%s_nx] ;
530 s_cxar[ii].i = 0.0f ;
531 }
532
533 csfft_cox( -1 , s_nfft , s_cxar ) ;
534
535 for( ii=1 ; ii < s_nftop ; ii++ ){
536 theta = TPI * (float)erand48(xran_pm) ;
537 CXROT(s_cxar[ii],theta) ;
538 s_cxar[s_nfft-ii].r = s_cxar[ii].r ;
539 s_cxar[s_nfft-ii].i = -s_cxar[ii].i ;
540 }
541 s_cxar[0].r = s_cxar[0].i ;
542 if( s_nfft2 > 0 ){
543 s_cxar[s_nfft2].i = 0.0f ;
544 if( nrand48(xran_pm)%2 == 1 )
545 s_cxar[s_nfft2].r = -s_cxar[s_nfft2].r ;
546 }
547
548 csfft_cox( 1 , s_nfft , s_cxar ) ;
549
550 for( ii=0 ; ii < s_nx ; ii++ ) ar[ii] = s_cxar[ii].r ;
551 THD_generic_retrend( s_nx , ar , -1 , 3 , s_pref , s_ffit ) ;
552 return ;
553 }
554
SORTS_FFTfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)555 static void SORTS_FFTfunc( double tzero, double tdelta ,
556 int npts, float ts[],
557 double ts_mean, double ts_slope,
558 void *ud, int nbriks, float *val )
559 {
560 int ii , nval ;
561
562 /** is this a "notification"? **/
563
564 if( val == NULL ){
565 #if 0
566 if( npts > 0 ){ /* the "start notification" */
567 } else { /* the "end notification" */
568 }
569 #endif
570 return ;
571 }
572
573 /** do the work **/
574
575 nval = MIN(nbriks,npts) ;
576 memcpy( val , ts , sizeof(float)*nval ) ;
577 float_scramble( nval , val ) ;
578 return ;
579 }
580