1 /*
2 afni/src/3dMSE.c
3 */
4
5 // Look for OpenMP macro
6 #ifdef USE_OMP
7 #include <omp.h>
8 #endif
9
10 // Include libraries
11 #include "mrilib.h"
12 #include <sys/mman.h>
13 #include <sys/types.h>
14
15 // Define constants
16 #define SPEARMAN 1
17 #define QUADRANT 2
18 #define PEARSON 3
19 #define ETA2 4
20
21 /* 3dMSE was created from 3dAutoTCorrelate by
22 Stan Colcombe and R. Cameron Craddock */
23
24 /*----------------------------------------------------------------*/
25 /**** Include these here for potential optimization for OpenMP ****/
26 /*----------------------------------------------------------------*/
27 /*! Pearson correlation of x[] and y[] (x and y are NOT modified.
28 And we know ahead of time that the time series have 0 mean
29 and L2 norm 1.
30 *//*--------------------------------------------------------------*/
31
zm_THD_pearson_corr(int n,float * x,float * y)32 float zm_THD_pearson_corr( int n, float *x , float *y ) /* inputs are */
33 { /* zero mean */
34 register float xy ; register int ii ; /* and norm=1 */
35 if( n%2 == 0 ){
36 xy = 0.0f ;
37 for( ii=0 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
38 } else {
39 xy = x[0]*y[0] ;
40 for( ii=1 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
41 }
42 return xy ;
43 }
44
45 /*----------------------------------------------------------------------------*/
vstep_print(void)46 static void vstep_print(void)
47 {
48 static int nn=0 ;
49 static char xx[10] = "0123456789" ;
50 fprintf(stderr , "%c" , xx[nn%10] ) ;
51 if( nn%10 == 9) fprintf(stderr,",") ;
52 nn++ ;
53 }
54
55 /*-----------------------------------------------------------------------------*/
sampleEn(float * xvec,long nval,long w,double r)56 double sampleEn( float* xvec, long nval, long w, double r )
57 {
58
59 long A = 0;
60 long B = 0;
61 double dist = 0.0;
62
63 long ii, jj, kk;
64
65 for ( ii=0; ii<(nval-w); ii++ )
66 {
67 dist = 0.0;
68 for ( jj=ii; jj<(nval-w); jj++ )
69 {
70 for ( kk=0; kk<w; kk++ )
71 {
72 /* lets make sure we aren't exceeding
73 arrays */
74 if ((ii+kk > nval) || (jj+kk > nval ))
75 {
76 WARNING_message("Indexing exceeds array bounds (%ld, %ld, %ld) (%s,%d)\n",
77 ii+kk,jj+kk,nval,__FILE__,__LINE__ );
78 }
79 else
80 {
81 dist = MAX( dist, fabs(xvec[ii+kk]-xvec[jj+kk]) );
82 }
83 }
84
85 if( dist < r )
86 {
87 B = B + 1;
88 }
89
90 if ((ii+kk+1 > nval) || (jj+kk+1 > nval ))
91 {
92 WARNING_message("Indexing exceeds array bounds (%ld, %ld, %ld) (%s,%d)\n",
93 ii+kk+1,jj+kk+1,nval,__FILE__,__LINE__ );
94 }
95 else if( MAX(dist, fabs(xvec[ii+kk+1]-xvec[jj+kk+1])) < r )
96 {
97 A = A + 1;
98 }
99 }
100 }
101
102 if (( A > 0 ) && ( B > 0 ))
103 {
104 dist = -1.0 * log((double)A / (double)B );
105 }
106 else
107 {
108 dist = 0.0;
109 }
110
111 return( dist );
112 }
113
114 /*----------------------------------------------------------------------------*/
115
main(int argc,char * argv[])116 int main( int argc , char *argv[] )
117 {
118 THD_3dim_dataset *xset , *cset, *mset=NULL ;
119 int nopt=1 , method=PEARSON , do_autoclip=0 ;
120 int nvox , nvals , ii, jj, kout, kin, polort=1 ;
121 int ix1,jy1,kz1, ix2, jy2, kz2 ;
122 char *prefix = "MSE" ;
123 byte *mask=NULL;
124 int nmask , abuc=1 ;
125 char str[32] , *cpt ;
126 int *imap = NULL ; MRI_vectim *xvectim ;
127
128 /* CC - number of scales to use in the calculation */
129 long num_scales = 5;
130
131 /* CC - the radius used in the SampleEn calculation */
132 long ent_win = 2;
133 double rthresh = 0.5;
134 double sd = 0.0;
135
136 /* CC - variables for writing out the results */
137 int nsubbriks = 0;
138 int subbrik = 0;
139 float * odset;
140
141 /* variables for holding results */
142 double * mse_results = NULL;
143
144 /*----*/
145
146 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
147
148 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
149 printf(
150 "Usage: 3dMSE [options] dset\n"
151 " Computes voxelwise multi-scale entropy."
152 "\n"
153 "Options:\n"
154 " -polort m = Remove polynomical trend of order 'm', for m=-1..3.\n"
155 " [default is m=1; removal is by least squares].\n"
156 " Using m=-1 means no detrending; this is only useful\n"
157 " for data/information that has been pre-processed.\n"
158 "\n"
159 " -autoclip = Clip off low-intensity regions in the dataset,\n"
160 " -automask = so that the correlation is only computed between\n"
161 " high-intensity (presumably brain) voxels. The\n"
162 " mask is determined the same way that 3dAutomask works.\n"
163 "\n"
164 " -mask mmm = Mask to define 'in-brain' voxels. Reducing the number\n"
165 " the number of voxels included in the calculation will\n"
166 " significantly speedup the calculation. Consider using\n"
167 " a mask to constrain the calculations to the grey matter\n"
168 " rather than the whole brain. This is also preferrable\n"
169 " to using -autoclip or -automask.\n"
170 "\n"
171 " -prefix p = Save output into dataset with prefix 'p', this file will\n"
172 " contain bricks for both 'weighted' or 'degree' centrality\n"
173 " [default prefix is 'MSE'].\n"
174 "\n"
175 " -scales N = The number of scales to be used in the calculation.\n"
176 " [default is 5].\n"
177 "\n"
178 " -entwin w = The window size used in the calculation.\n"
179 " [default is 2].\n"
180 "\n"
181 " -rthresh r = The radius threshold for determining if values are the\n"
182 " same in the SampleEn calculation, in fractions of the\n"
183 " standard deviation.\n"
184 " [default is .5].\n"
185 "\n"
186 "Notes:\n"
187 " * The output dataset is a bucket type of floats.\n"
188 "\n"
189 "-- RWCox - 31 Jan 2002 and 16 Jul 2010\n"
190 "-- Cameron Craddock - 26 Sept 2015 \n"
191 ) ;
192 PRINT_AFNI_OMP_USAGE("3dMSE",NULL) ;
193 PRINT_COMPILE_DATE ; exit(0) ;
194 }
195
196 mainENTRY("3dMSE main"); machdep(); PRINT_VERSION("3dMSE");
197 AFNI_logger("3dMSE",argc,argv);
198
199
200 /*-- option processing --*/
201
202 while( nopt < argc && argv[nopt][0] == '-' ){
203
204 if( strcmp(argv[nopt],"-time") == 0 ){
205 abuc = 0 ; nopt++ ; continue ;
206 }
207
208 if( strcmp(argv[nopt],"-autoclip") == 0 ||
209 strcmp(argv[nopt],"-automask") == 0 ){
210
211 do_autoclip = 1 ; nopt++ ; continue ;
212 }
213
214 if( strcmp(argv[nopt],"-mask") == 0 ){
215 mset = THD_open_dataset(argv[++nopt]);
216 CHECK_OPEN_ERROR(mset,argv[nopt]);
217 nopt++ ; continue ;
218 }
219
220 if( strcmp(argv[nopt],"-prefix") == 0 ){
221 prefix = strdup(argv[++nopt]) ;
222 if( !THD_filename_ok(prefix) ){
223 ERROR_exit("Illegal value after -prefix!") ;
224 }
225 nopt++ ; continue ;
226 }
227
228 if( strcmp(argv[nopt],"-scales") == 0 ){
229 long val = (long)strtod(argv[++nopt],&cpt) ;
230 if( *cpt != '\0' || val < 0.0 ){
231 ERROR_exit("Illegal value (%ld) after -scales!", val) ;
232 }
233 num_scales = val ; nopt++ ; continue ;
234 }
235
236 if( strcmp(argv[nopt],"-rthresh") == 0 ){
237 double val = (double)strtod(argv[++nopt],&cpt) ;
238 if( *cpt != '\0' || val < 0.0 ){
239 ERROR_exit("Illegal value (%f) after -rthresh!", val) ;
240 }
241 rthresh = val ; nopt++ ; continue ;
242 }
243
244 if( strcmp(argv[nopt],"-entwin") == 0 ){
245 int val = (int)strtod(argv[++nopt],&cpt) ;
246 if( *cpt != '\0' || val < 0 ){
247 ERROR_exit("Illegal value after -entwin!") ;
248 }
249 ent_win = val ; nopt++ ; continue ;
250 }
251
252 if( strcmp(argv[nopt],"-polort") == 0 ){
253 int val = (int)strtod(argv[++nopt],&cpt) ;
254 if( *cpt != '\0' || val < -1 || val > 3 ){
255 ERROR_exit("Illegal value after -polort!") ;
256 }
257 polort = val ; nopt++ ; continue ;
258 }
259
260 ERROR_exit("Illegal option: %s",argv[nopt]) ;
261 }
262
263 /*-- open dataset, check for legality --*/
264
265 if( nopt >= argc ) ERROR_exit("Need a dataset on command line!?") ;
266
267 xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);
268
269 if( DSET_NVALS(xset) < ent_win * num_scales + 1 )
270 ERROR_exit("Input dataset %s does not have enough sub-bricks!",argv[nopt]) ;
271 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
272
273 /*-- compute mask array, if desired --*/
274 nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
275
276 /* if a mask was specified make sure it is appropriate */
277 if( mset ){
278
279 if( DSET_NVOX(mset) != nvox )
280 ERROR_exit("Input and mask dataset differ in number of voxels!") ;
281
282 /* make a mask */
283 mask = THD_makemask(mset, 0, 1.0, 0.0) ;
284
285 /* determine number of voxels in the mask */
286 nmask = THD_countmask( nvox , mask ) ;
287
288 INFO_message("%d voxels in -mask dataset",nmask) ;
289 if( nmask < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmask);
290
291 /* update running memory statistics to reflect loading the image */
292 DSET_unload(mset) ;
293 }
294 /* if automasking is requested, handle that now */
295 else if( do_autoclip ){
296 mask = THD_automask( xset ) ;
297 nmask = THD_countmask( nvox , mask ) ;
298 INFO_message("%d voxels survive -autoclip",nmask) ;
299 if( nmask < 2 ) ERROR_exit("Only %d voxels in -automask!",nmask);
300 }
301 /* otherwise we use all of the voxels in the image */
302 else {
303 nmask = nvox ;
304 INFO_message("computing for all %d voxels",nmask) ;
305 }
306
307 /*-- create vectim from input dataset --*/
308 INFO_message("vectim-izing input dataset") ;
309
310 /*-- CC added in mask to reduce the size of xvectim -- */
311 xvectim = THD_dset_to_vectim( xset , mask , 0 ) ;
312 if( xvectim == NULL ) ERROR_exit("Can't create vectim?!") ;
313
314 /*--- CC the vectim contains a mapping between voxel index and mask index,
315 tap into that here to avoid duplicating memory usage ---*/
316
317 if( mask != NULL )
318 {
319 imap = xvectim->ivec;
320
321 /* --- CC free the mask */
322 free(mask); mask=NULL;
323 }
324
325 /* -- CC unloading the dataset to reduce memory usage ?? -- */
326 DSET_unload(xset) ;
327
328 /* detrend the data */
329 if( polort >= 0 )
330 {
331 for( ii=0 ; ii < xvectim->nvec ; ii++ )
332 { /* remove polynomial trend */
333 DETREND_polort(polort,nvals,VECTIM_PTR(xvectim,ii)) ;
334 }
335 }
336
337 /* -- this procedure does not change time series that have zero variance -- */
338 THD_vectim_normalize(xvectim) ; /* L2 norm = 1 */
339
340 /* after this normalization, the sd = sqrt(1/nvals) */
341 sd = sqrt(1.0/(double)xvectim->nvals);
342
343 /* -- allocate memory to hold the results as they are being calculated -- */
344 if( ( mse_results = (double*)calloc( num_scales*nmask, sizeof(double) )) == NULL )
345 {
346 ERROR_message( "Could not allocate %d byte array for MSE calculation\n",
347 nmask*sizeof(long));
348 }
349
350 /*-- tell the user what we are about to do --*/
351 INFO_message( "Calculating multi-scale entropy for %d scales, window length = %ld, and threshold r = %f.\n",
352 num_scales,ent_win,rthresh);
353
354 /*---------- loop over mask voxels, correlate ----------*/
355 AFNI_OMP_START ;
356 #pragma omp parallel if( nmask > 999 )
357 {
358 long ithr, nthr;
359 long lout, lin, kin, scale;
360 float *xsar;
361 float *ysar = NULL;
362 float *ysar_orig = NULL;
363 long vstep, vii;
364 long scale_nvals;
365
366 /*-- get information about who we are --*/
367 #ifdef USE_OMP
368 ithr = omp_get_thread_num() ;
369 nthr = omp_get_num_threads() ;
370 if( ithr == 0 ) INFO_message("%d OpenMP threads started",nthr) ;
371 #else
372 ithr = 0 ; nthr = 1 ;
373 #endif
374
375 /*-- For the progress tracker, we want to print out 50 numbers,
376 figure out a number of loop iterations that will make this easy */
377 vstep = (long)( nmask / (nthr*50.0f) + 0.901f ) ; vii = 0 ;
378 if(ithr == 0 ) fprintf(stderr,"Looping:") ;
379
380 /* allocate memory */
381 #pragma omp critical (MALLOC)
382 {
383 ysar = (float*)malloc( xvectim->nvals * sizeof(float));
384 }
385
386 if( ysar == NULL )
387 {
388 WARNING_message("Could not allocate scratch memory for thread %d (%s,%d)\n",
389 ithr,__FILE__,__LINE__ );
390 }
391
392 ysar_orig = ysar;
393 #pragma omp for schedule(static, 1)
394 for( lout=0 ; lout < xvectim->nvec ; lout++ ) /*----- outer voxel loop -----*/
395 {
396 if( ysar != NULL )
397 {
398 if( ithr == 0 && vstep > 2 ) /* allow small dsets 16 Jun 2011 [rickr] */
399 {
400 vii++;
401 if( vii%vstep == vstep/2 ) vstep_print();
402 }
403
404 /* get ref time series from this voxel */
405 xsar = VECTIM_PTR(xvectim,lout) ;
406
407 /* calculate the entropy at scale = 1 */
408 mse_results[ lout ] =
409 sampleEn( xsar, xvectim->nvals, ent_win, rthresh * sd );
410
411 /* calculate the entropy at the remaining scales */
412 for( scale = 2; scale <= num_scales; scale++ )
413 {
414 scale_nvals = (long)floor((double)xvectim->nvals / (double)scale);
415 for(lin=0; lin < scale_nvals; lin++ )
416 {
417 if( lin > xvectim->nvals )
418 {
419 WARNING_message("Avoiding ysar overflow %d > %d (%s,%d)\n",
420 lin,xvectim->nvals,__FILE__,__LINE__ );
421 }
422 else
423 {
424 ysar[ lin ] = 0.0;
425 for(kin=0; kin<scale; kin++ )
426 {
427 ysar[ lin ] = ysar[ lin ] + xsar[ lin*scale + kin ];
428 }
429 ysar[ lin ] = ysar[ lin ] / (double)scale;
430 }
431 }
432 mse_results[ ((scale-1)*xvectim->nvec) + lout ] =
433 sampleEn( ysar, scale_nvals, ent_win, rthresh * sd );
434 }
435 }
436 } /* end of outer loop over ref voxels */
437
438 if( ithr == 0 ) fprintf(stderr,".\n") ;
439
440 #pragma omp critical (MALLOC)
441 {
442 if( ysar != NULL )
443 {
444 free(ysar);
445 ysar = NULL;
446 }
447 }
448
449 } /* end OpenMP */
450 AFNI_OMP_END ;
451
452 /* update the user so that they know what we are up to */
453 #ifdef USE_OMP
454 INFO_message ("AFNI_OMP finished\n");
455 #endif
456
457 /*---------- Finish up ---------*/
458
459 /*-- create output dataset --*/
460 cset = EDIT_empty_copy( xset ) ;
461
462 /*-- configure the output dataset */
463 if( abuc ){
464 EDIT_dset_items( cset ,
465 ADN_prefix , prefix ,
466 ADN_nvals , num_scales , /* subbricks */
467 ADN_ntt , 0 , /* no time axis */
468 ADN_type , HEAD_ANAT_TYPE ,
469 ADN_func_type , ANAT_BUCK_TYPE ,
470 ADN_datum_all , MRI_float ,
471 ADN_none ) ;
472 } else {
473 EDIT_dset_items( cset ,
474 ADN_prefix , prefix ,
475 ADN_nvals , num_scales , /* subbricks */
476 ADN_ntt , num_scales , /* num times */
477 ADN_ttdel , 1.0 , /* fake TR */
478 ADN_nsl , 0 , /* no slice offsets */
479 ADN_type , HEAD_ANAT_TYPE ,
480 ADN_func_type , ANAT_EPI_TYPE ,
481 ADN_datum_all , MRI_float ,
482 ADN_none ) ;
483 }
484
485 /* add history information to the hearder */
486 tross_Make_History( "3dMSE" , argc,argv , cset ) ;
487
488 ININFO_message("Creating output dataset in memory") ;
489
490 /* -- Configure the subbriks and copy in the data */
491 for (subbrik=0; subbrik<num_scales; subbrik++ )
492 {
493 EDIT_BRICK_TO_NOSTAT(cset,subbrik) ; /* stat params */
494 /* CC this sets the subbrik scaling factor, which we will probably want
495 to do again after we calculate the voxel values */
496 EDIT_BRICK_FACTOR(cset,subbrik,1.0) ; /* scale factor */
497
498 sprintf(str,"Multi-scale entropy (%d)", subbrik) ;
499
500 EDIT_BRICK_LABEL(cset,subbrik,str) ;
501 EDIT_substitute_brick(cset,subbrik,MRI_float,NULL) ; /* make array */
502
503 /* copy measure data into the subbrik */
504 odset = DSET_ARRAY(cset,subbrik);
505
506 for( kout = 0; kout < nmask; kout++ )
507 {
508 if ( imap != NULL )
509 {
510 ii = imap[kout] ; /* ii= source voxel (we know that ii is in the mask) */
511 }
512 else
513 {
514 ii = kout ;
515 }
516
517 if( mse_results == NULL )
518 {
519 WARNING_message("MSE Results is NULL %d > %d (%s,%d)\n",
520 ii,DSET_NVOX(cset),__FILE__,__LINE__ );
521 }
522 else if( ii >= DSET_NVOX(cset) )
523 {
524 WARNING_message("Avoiding odset overflow %d > %d (%s,%d)\n",
525 ii,DSET_NVOX(cset),__FILE__,__LINE__ );
526 }
527 else if( (subbrik*nmask)+kout >= num_scales*nmask )
528 {
529 WARNING_message("Avoiding mse_results overflow %d >= %d (%s,%d)\n",
530 (subbrik*nmask)+kout,num_scales*nmask,__FILE__,__LINE__ );
531 }
532 else
533 {
534 odset[ ii ] = (double)(mse_results[(subbrik*nmask)+kout]);
535 }
536 }
537 }
538
539 /* we are done with this memory, and can kill it now*/
540 if(mse_results)
541 {
542 free(mse_results);
543 mse_results=NULL;
544 }
545
546 INFO_message("Done..\n") ;
547
548 /* toss some trash */
549 VECTIM_destroy(xvectim) ;
550 DSET_delete(xset) ;
551
552 /* finito */
553 INFO_message("Writing output dataset to disk [%s bytes]",
554 commaized_integer_string(cset->dblk->total_bytes)) ;
555
556 /* write the dataset */
557 DSET_write(cset) ;
558 WROTE_DSET(cset) ;
559
560 exit(0) ;
561 }
562