1 /*
2 afni/src/3dLFCD.c
3 */
4
5 /* 3dLFCD was created from 3dAutoTCorrelate by
6 R. Cameron Craddock */
7
8 // Look for OpenMP macro
9 #ifdef USE_OMP
10 #include <omp.h>
11 #endif
12
13 // Include libraries
14 #include "mrilib.h"
15 #include <sys/mman.h>
16 #include <sys/types.h>
17
18 // Define constants
19 #define SPEARMAN 1
20 #define QUADRANT 2
21 #define PEARSON 3
22 #define ETA2 4
23
24 #define MAX_NUM_TAGS 32
25 #define MAX_TAG_LEN 256
26
27 #define FACES 0
28 #define FACES_EDGES 1
29 #define FACES_EDGES_CORNERS 2
30
31 #define X 0
32 #define Y 1
33 #define Z 2
34
35 #define faces_neighborhood_num 6
36 int faces_neighborhood[18] = {
37 0, 0, 1,
38 0, 0, -1,
39 0, 1, 0,
40 0, -1, 0,
41 1, 0, 0,
42 -1, 0, 0
43 };
44
45 /* x y z */
46 #define faces_edges_neighborhood_num 18
47 const int faces_edges_neighborhood[54] = {
48 0, -1, -1,
49 -1, 0, -1,
50 0, 0, -1,
51 1, 0, -1,
52 0, 1, -1,
53 -1, -1, 0,
54 -1, 0, 0,
55 -1, 1, 0,
56 0, -1, 0,
57 0, 1, 0,
58 1, -1, 0,
59 1, 0, 0,
60 1, 1, 0,
61 0, -1, 1,
62 -1, 0, 1,
63 0, 0, 1,
64 1, 0, 1,
65 0, 1, 1
66 };
67
68 /* x y z */
69 #define faces_edges_corners_neighborhood_num 26
70 const int faces_edges_corners_neighborhood[78] = {
71 -1, -1, -1,
72 -1, 0, -1,
73 -1, 1, -1,
74 0, -1, -1,
75 0, 0, -1,
76 0, 1, -1,
77 1, -1, -1,
78 1, 0, -1,
79 1, 1, -1,
80 -1, -1, 0,
81 -1, 0, 0,
82 -1, 1, 0,
83 0, -1, 0,
84 0, 1, 0,
85 1, -1, 0,
86 1, 0, 0,
87 1, 1, 0,
88 -1, -1, 1,
89 -1, 0, 1,
90 -1, 1, 1,
91 0, -1, 1,
92 0, 0, 1,
93 0, 1, 1,
94 1, -1, 1,
95 1, 0, 1,
96 1, 1, 1,
97 };
98
99 /* CC macro for updating mem stats */
100 #define INC_MEM_STATS( INC, TAG ) \
101 { \
102 if( MEM_PROF == 1 ) \
103 { \
104 int ndx = 0; \
105 while( ndx < mem_num_tags ) \
106 { \
107 if( strncmp( mem_tags[ndx], TAG, MAX_TAG_LEN ) == 0 ) \
108 { \
109 break; \
110 } \
111 ndx++; \
112 } \
113 if(( ndx >= mem_num_tags ) && (ndx < MAX_NUM_TAGS)) \
114 { \
115 /* adding a new tag */ \
116 strncpy( mem_tags[ ndx ], TAG, (MAX_TAG_LEN-1) ); \
117 mem_allocated[ ndx ] = 0; \
118 mem_freed[ ndx ] = 0; \
119 mem_num_tags++; \
120 } \
121 if( ndx < MAX_NUM_TAGS ) \
122 { \
123 mem_allocated[ ndx ] += (long)(INC); \
124 if ((long)(INC) > 1024 ) WARNING_message("Incrementing memory for %s by %ldB\n", TAG, (INC)); \
125 } \
126 else WARNING_message("No room in mem profiler for %s\n", TAG ); \
127 } \
128 total_mem += (long)(INC); \
129 running_mem += (long)(INC); \
130 if (running_mem > peak_mem) peak_mem = running_mem; \
131 }
132
133 #define DEC_MEM_STATS( DEC, TAG ) \
134 { \
135 if( MEM_PROF == 1 ) \
136 { \
137 int ndx = 0; \
138 while( ndx < mem_num_tags ) \
139 { \
140 if( strncmp( mem_tags[ndx], TAG, MAX_TAG_LEN ) == 0 ) \
141 { \
142 break; \
143 } \
144 else ndx++ ; \
145 } \
146 if(( ndx >= mem_num_tags ) && (ndx < MAX_NUM_TAGS)) \
147 { \
148 WARNING_message("Could not find tag %s in mem profiler\n", TAG ); \
149 } \
150 else \
151 { \
152 mem_freed[ ndx ] += (long)(DEC); \
153 if ((long)(DEC) > 1024 ) INFO_message("Free %ldB of memory for %s\n", (DEC), TAG); \
154 } \
155 } \
156 running_mem -= (long)(DEC); \
157 }
158
159 #define PRINT_MEM_STATS( TAG ) \
160 if ( MEM_STAT == 1 ) \
161 { \
162 INFO_message("\n======\n== Mem Stats (%s): Running %fB, Total %3.3fMB, Peak %3.3fMB\n", \
163 TAG, \
164 (double)(running_mem), \
165 (double)(total_mem/(1024.0*1024.0)), \
166 (double)(peak_mem/(1024.0*1024.0))); \
167 INFO_message("\n======\n== Mem Stats (%s): Running %3.3fMB, Total %3.3fMB, Peak %3.3fMB\n", \
168 TAG, \
169 (double)(running_mem/(1024.0*1024.0)), \
170 (double)(total_mem/(1024.0*1024.0)), \
171 (double)(peak_mem/(1024.0*1024.0))); \
172 if( MEM_PROF == 1 ) \
173 { \
174 int ndx = 0; \
175 INFO_message("== Memory Profile\n"); \
176 for( ndx=0; ndx < mem_num_tags; ndx++ ) \
177 { \
178 INFO_message("%s: %ld allocated %ld freed\n", mem_tags[ndx], \
179 mem_allocated[ndx], mem_freed[ndx] ); \
180 } \
181 } \
182 }
183
184
185 /* CC define nodes for list that will be used
186 for region growing */
187 typedef struct _list_node list_node;
188
189 // Define list node structure
190 struct _list_node
191 {
192 long vox_vol_ndx;
193 long ix;
194 long jy;
195 long kz;
196
197 list_node* next;
198 };
199
200
201 /* free the list of list_nodes */
202 #define FREE_LIST( list ) \
203 { \
204 list_node *pptr; \
205 list_node *hptr = list; \
206 \
207 while(hptr != NULL ) \
208 { \
209 /* increment node pointers */ \
210 pptr = hptr; \
211 hptr = hptr->next; \
212 \
213 /* delete the node */ \
214 if(pptr != NULL) \
215 { \
216 /* -- update running memory estimate to reflect memory allocation */ \
217 DEC_MEM_STATS(( sizeof(list_node)), "list nodes"); \
218 free(pptr); \
219 pptr=NULL; \
220 } \
221 } \
222 list = NULL; \
223 }
224
225 /* freeing all of the allocated mem on an error can get a little messy. instead
226 we can use this macro to check what has been allocated and kill it. this of
227 course requires strict discipline for initiazing all pointers to NULL and
228 resetting them to NULL when free'd. i should be able to handle that */
229 #define CHECK_AND_FREE_ALL_ALLOCATED_MEM \
230 { \
231 /* eliminate DSETS */ \
232 if ( mset != NULL ) \
233 { \
234 DSET_unload(mset) ; \
235 DSET_delete(mset) ; \
236 mset = NULL ; \
237 } \
238 \
239 if ( xset != NULL ) \
240 { \
241 DSET_unload(xset) ; \
242 DSET_delete(xset) ; \
243 xset = NULL ; \
244 } \
245 \
246 if ( cset != NULL ) \
247 { \
248 DSET_unload(cset) ; \
249 DSET_delete(cset) ; \
250 cset = NULL ; \
251 } \
252 \
253 /* free the xvectim */ \
254 if ( xvectim != NULL ) \
255 { \
256 VECTIM_destroy(xvectim) ; \
257 xvectim = NULL ; \
258 } \
259 \
260 /* free allocated mems */ \
261 if( mask != NULL ) \
262 { \
263 free(mask) ; \
264 mask = NULL ; \
265 } \
266 \
267 if( vol_ndx_to_mask_ndx != NULL ) \
268 { \
269 free(vol_ndx_to_mask_ndx) ; \
270 vol_ndx_to_mask_ndx = NULL ; \
271 } \
272 }
273
274 /* being good and cleaning up before erroring out can be a pain, and messy, lets
275 make a variadic macro that does it for us. */
276 #define ERROR_EXIT_CC( ... ) \
277 { \
278 CHECK_AND_FREE_ALL_ALLOCATED_MEM; \
279 ERROR_exit( __VA_ARGS__ ); \
280 }
281
282
283
284
285 /*----------------------------------------------------------------*/
286 /**** Include these hesre for potential optimization for OpenMP ****/
287 /*----------------------------------------------------------------*/
288 /*! Pearson correlation of x[] and y[] (x and y are NOT modified.
289 And we know ahead of time that the time series have 0 mean
290 and L2 norm 1.
291 *//*--------------------------------------------------------------*/
292
zm_THD_pearson_corr(int n,float * x,float * y)293 float zm_THD_pearson_corr( int n, float *x , float *y ) /* inputs are */
294 { /* zero mean */
295 register float xy ; register int ii ; /* and norm=1 */
296 if( n%2 == 0 ){
297 xy = 0.0f ;
298 for( ii=0 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
299 } else {
300 xy = x[0]*y[0] ;
301 for( ii=1 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
302 }
303 return xy ;
304 }
305
306
307 /* */
308
309 /*----------------------------------------------------------------*/
310 /* General correlation calculation. */
311
312 #if 0
313 float my_THD_pearson_corr( int n, float *x , float *y )
314 {
315 float xv,yv,xy , vv,ww , xm,ym ;
316 register int ii ;
317
318 xm = ym = 0.0f ;
319 for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
320 xm /= n ; ym /= n ;
321 xv = yv = xy = 0.0f ;
322 for( ii=0 ; ii < n ; ii++ ){
323 vv = x[ii]-xm ; ww = y[ii]-ym ; xv += vv*vv ; yv += ww*ww ; xy += vv*ww ;
324 }
325
326 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
327 return xy/sqrtf(xv*yv) ;
328 }
329 #endif
330
331 /*----------------------------------------------------------------*/
332 /*! eta^2 (Cohen, NeuroImage 2008) 25 Jun 2010 [rickr]
333 *
334 * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ]
335 * ------------------------------------
336 * SUM[ (a_i - M )^2 + (b_i - M )^2 ]
337 *
338 * where o a_i and b_i are the vector elements
339 * o m_i = (a_i + b_i)/2
340 * o M = mean across both vectors
341 -----------------------------------------------------------------*/
342
my_THD_eta_squared(int n,float * x,float * y)343 float my_THD_eta_squared( int n, float *x , float *y )
344 {
345 float num , denom , gm , lm, vv, ww;
346 register int ii ;
347
348 gm = 0.0f ;
349 for( ii=0 ; ii < n ; ii++ ){ gm += x[ii] + y[ii] ; }
350 gm /= (2.0f*n) ;
351
352 num = denom = 0.0f ;
353 for( ii=0 ; ii < n ; ii++ ){
354 lm = 0.5f * ( x[ii] + y[ii] ) ;
355 vv = (x[ii]-lm); ww = (y[ii]-lm); num += ( vv*vv + ww*ww );
356 vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww );
357 }
358
359 if( num < 0.0f || denom <= 0.0f || num >= denom ) return 0.0f ;
360 return (1.0f - num/denom) ;
361 }
362
363 /*----------------------------------------------------------------------------*/
vstep_print(void)364 static void vstep_print(void)
365 {
366 static int nn=0 ;
367 static char xx[10] = "0123456789" ;
368 fprintf(stderr , "%c" , xx[nn%10] ) ;
369 if( nn%10 == 9) fprintf(stderr,",") ;
370 nn++ ;
371 }
372
373 /*-----------------------------------------------------------------------------*/
374
375 /*----------------------------------------------------------------------------*/
376
main(int argc,char * argv[])377 int main( int argc , char *argv[] )
378 {
379 THD_3dim_dataset *xset = NULL;
380 THD_3dim_dataset *cset = NULL;
381 THD_3dim_dataset *mset = NULL ;
382 int nopt=1 , method=PEARSON , do_autoclip=0 ;
383 int nvox , nvals , ii, polort=1 ;
384 char *prefix = "LFCD" ;
385 byte *mask=NULL;
386 int nmask , abuc=1 ;
387 char str[32] , *cpt = NULL;
388 int *mask_ndx_to_vol_ndx = NULL;
389 MRI_vectim *xvectim = NULL ;
390 float (*corfun)(int,float *,float*) = NULL ;
391
392 /* CC - we will have two subbricks: binary and weighted centrality */
393 int nsubbriks = 2;
394 int subbrik = 0;
395 float * bodset;
396 float * wodset;
397
398 /* CC - added flags for thresholding correlations */
399 double thresh = 0.0;
400
401 double eps = 1e-9;
402
403 /* CC - defines the neighborhood */
404 int num_neighbors = faces_neighborhood_num;
405 const int * neighborhood = faces_neighborhood;
406
407 /* CC - variables to assist going back and forth between mask and volume indices */
408 long * vol_ndx_to_mask_ndx = NULL;
409
410 /* CC - variables for tracking memory usage stats */
411 char mem_tags[MAX_NUM_TAGS][MAX_TAG_LEN];
412 long mem_allocated[MAX_NUM_TAGS];
413 long mem_freed[MAX_NUM_TAGS];
414 long mem_num_tags = 0;
415 long running_mem = 0;
416 long peak_mem = 0;
417 long total_mem = 0;
418 int MEM_PROF = 0;
419 int MEM_STAT = 0;
420 /*----*/
421
422 /* the number of zero variance voxels removed from the data */
423 int rem_count = 0;
424
425 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
426
427 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
428 printf(
429 "Usage: 3dLFCD [options] dset\n"
430 " Computes voxelwise local functional connectivity density as defined in:\n"
431 " Tomasi, D and Volkow, PNAS, May 2010, 107 (21) 9885-9890; \n"
432 " DOI: 10.1073/pnas.1001414107 \n\n"
433 " The results are stored in a new 3D bucket dataset\n as floats to preserve\n"
434 " their values. Local functional connectivity density (LFCD; as opposed to global\n"
435 " functional connectivity density, see 3dDegreeCentrality), reflects\n"
436 " the extent of the correlation of a voxel within its locally connected cluster.\n\n"
437 " Conceptually the process involves: \n"
438 " 1. Calculating the correlation between voxel time series for\n"
439 " every pair of voxels in the brain (as determined by masking)\n"
440 " 2. Applying a threshold to the resulting correlations to exclude\n"
441 " those that might have arisen by chance\n"
442 " 3. Find the cluster of above-threshold voxels that are spatially\n"
443 " connected to the target voxel.\n"
444 " 4. Count the number of voxels in the local cluster.\n"
445 " Practically the algorithm is ordered differently to optimize for\n"
446 " computational time and memory usage.\n\n"
447 " The procedure described in the paper defines a voxels\n"
448 " neighborhood to be the 6 voxels with which it shares a face.\n"
449 " This definition can be changed to include edge and corner \n"
450 " voxels using the -neighborhood flags below.\n\n"
451 " LFCD is a localized variant of binarized degree centrality,\n"
452 " the weighted alternative is calculated by changing step 4\n"
453 " above to calculate the sum of the correlation coefficients\n"
454 " between the seed region and the neigbors. 3dLFCD outputs\n"
455 " both of these values (in seperate briks), since they are\n"
456 " so easy to calculate in tandem.\n\n"
457 " You might prefer to calculate this on your data after\n"
458 " spatial normalization, so that the range of values are\n"
459 " consistent between datatsets. Similarly the same brain mask\n"
460 " should be used for all datasets that will be directly compared.\n\n"
461 " The original paper used a correlation threshold = 0.6 and \n"
462 " excluded all voxels with tSNR < 50. 3dLFCD does not discard\n"
463 " voxels based on tSNR, this would need to be done beforehand.\n"
464 "\n"
465 "Options:\n"
466 " -pearson = Correlation is the normal Pearson (product moment)\n"
467 " correlation coefficient [default].\n"
468 #if 0
469 " -spearman = Correlation is the Spearman (rank) correlation\n"
470 " coefficient.\n"
471 " -quadrant = Correlation is the quadrant correlation coefficient.\n"
472 #else
473 " -spearman AND -quadrant are disabled at this time :-(\n"
474 #endif
475 "\n"
476 " -thresh r = exclude correlations <= r from calculations\n"
477 "\n"
478 " -faces = define neighborhood to include face touching\n"
479 " edges (default)\n"
480 " -faces_edges = define neighborhood to include face and \n"
481 " edge touching voxels\n"
482 " -faces_edges_corners = define neighborhood to include face, \n"
483 " edge, and corner touching voxels\n"
484 "\n"
485 " -polort m = Remove polynomial trend of order 'm', for m=-1..3.\n"
486 " [default is m=1; removal is by least squares].\n"
487 " Using m=-1 means no detrending; this is only useful\n"
488 " for data/information that has been pre-processed.\n"
489 "\n"
490 " -autoclip = Clip off low-intensity regions in the dataset,\n"
491 " -automask = so that the correlation is only computed between\n"
492 " high-intensity (presumably brain) voxels. The\n"
493 " mask is determined the same way that 3dAutomask works.\n"
494 " This is done automatically if no mask is proveded.\n"
495 "\n"
496 " -mask mmm = Mask to define 'in-brain' voxels. Reducing the number\n"
497 " the number of voxels included in the calculation will\n"
498 " significantly speedup the calculation. Consider using\n"
499 " a mask to constrain the calculations to the grey matter\n"
500 " rather than the whole brain. This is also preferrable\n"
501 " to using -autoclip or -automask.\n"
502 "\n"
503 " -prefix p = Save output into dataset with prefix 'p', this file will\n"
504 " contain bricks for both 'weighted' and 'binarized' lFCD\n"
505 " [default prefix is 'LFCD'].\n"
506 "\n"
507 "Notes:\n"
508 " * The output dataset is a bucket type of floats.\n"
509 " * The program prints out an estimate of its memory used\n"
510 " when it ends. It also prints out a progress 'meter'\n"
511 " to keep you pacified.\n"
512 "\n"
513 "-- RWCox - 31 Jan 2002 and 16 Jul 2010\n"
514 "-- Cameron Craddock - 13 Nov 2015 \n"
515 ) ;
516 PRINT_AFNI_OMP_USAGE("3dLFCD",NULL) ;
517 PRINT_COMPILE_DATE ; exit(0) ;
518 }
519
520 mainENTRY("3dLFCD main"); machdep(); PRINT_VERSION("3dLFCD");
521 AFNI_logger("3dLFCD",argc,argv);
522
523 /*-- option processing --*/
524
525 while( nopt < argc && argv[nopt][0] == '-' ){
526
527 if( strcmp(argv[nopt],"-time") == 0 ){
528 abuc = 0 ; nopt++ ; continue ;
529 }
530
531 if( strcmp(argv[nopt],"-autoclip") == 0 ||
532 strcmp(argv[nopt],"-automask") == 0 ){
533
534 do_autoclip = 1 ; nopt++ ; continue ;
535 }
536
537 if( strcmp(argv[nopt],"-mask") == 0 ){
538 /* mset is opened here, but not loaded? */
539 mset = THD_open_dataset(argv[++nopt]);
540 CHECK_OPEN_ERROR(mset,argv[nopt]);
541 nopt++ ; continue ;
542 }
543
544 if( strcmp(argv[nopt],"-pearson") == 0 ){
545 method = PEARSON ; nopt++ ; continue ;
546 }
547
548 #if 0
549 if( strcmp(argv[nopt],"-spearman") == 0 ){
550 method = SPEARMAN ; nopt++ ; continue ;
551 }
552
553 if( strcmp(argv[nopt],"-quadrant") == 0 ){
554 method = QUADRANT ; nopt++ ; continue ;
555 }
556 #endif
557
558 if( strcmp(argv[nopt],"-eta2") == 0 ){
559 method = ETA2 ; nopt++ ; continue ;
560 }
561
562 if( strcmp(argv[nopt],"-prefix") == 0 ){
563 prefix = strdup(argv[++nopt]) ;
564 if( !THD_filename_ok(prefix) ){
565 ERROR_EXIT_CC("Illegal value after -prefix!") ;
566 }
567 nopt++ ; continue ;
568 }
569
570 if( strcmp(argv[nopt],"-thresh") == 0 ){
571 double val = (double)strtod(argv[++nopt],&cpt) ;
572 if( *cpt != '\0' || val >= 1.0 || val < 0.0 ){
573 ERROR_EXIT_CC("Illegal value (%f) after -thresh!", val) ;
574 }
575 thresh = val ; nopt++ ; continue ;
576 }
577
578 if( strcmp(argv[nopt],"-faces") == 0 ){
579 num_neighbors = faces_neighborhood_num;
580 neighborhood = faces_neighborhood;
581 nopt++ ; continue ;
582 }
583
584 if( strcmp(argv[nopt],"-faces_edges") == 0 ){
585 num_neighbors = faces_edges_neighborhood_num;
586 neighborhood = faces_edges_neighborhood;
587 nopt++ ; continue ;
588 }
589
590 if( strcmp(argv[nopt],"-faces_edges_neighborhood") == 0 ){
591 num_neighbors = faces_edges_corners_neighborhood_num;
592 neighborhood = faces_edges_corners_neighborhood;
593 nopt++ ; continue ;
594 }
595
596 if( strcmp(argv[nopt],"-polort") == 0 ){
597 int val = (int)strtod(argv[++nopt],&cpt) ;
598 if( *cpt != '\0' || val < -1 || val > 3 ){
599 ERROR_EXIT_CC("Illegal value after -polort!") ;
600 }
601 polort = val ; nopt++ ; continue ;
602 }
603
604 if( strcmp(argv[nopt],"-mem_stat") == 0 ){
605 MEM_STAT = 1 ; nopt++ ; continue ;
606 }
607
608 if( strncmp(argv[nopt],"-mem_profile",8) == 0 ){
609 MEM_PROF = 1 ; nopt++ ; continue ;
610 }
611
612 ERROR_EXIT_CC("Illegal option: %s",argv[nopt]) ;
613 }
614
615 /*-- open dataset, check for legality --*/
616
617 if( nopt >= argc ) ERROR_EXIT_CC("Need a dataset on command line!?") ;
618
619 INFO_message("Calculating LFCD from %s using mask, r_threshold %lf and %d voxels in the neighborhood\n",
620 argv[nopt], thresh, num_neighbors);
621
622 xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);
623
624 if( DSET_NVALS(xset) < 3 )
625 ERROR_EXIT_CC("Input dataset %s does not have 3 or more sub-bricks!",
626 argv[nopt]) ;
627 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
628
629 /*-- compute mask array, if desired --*/
630 nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
631 INC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");
632 PRINT_MEM_STATS("inset");
633
634 /* if a mask was specified make sure it is appropriate */
635 if( mset ){
636
637 if( DSET_NVOX(mset) != nvox )
638 ERROR_EXIT_CC("Input and mask dataset differ in number of voxels!") ;
639 mask = THD_makemask(mset, 0, 1.0, 0.0) ;
640
641 /* update running memory statistics to reflect loading the image */
642 INC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );
643 PRINT_MEM_STATS( "mset load" );
644
645 /* iupdate statistics to reflect creating mask array */
646 nmask = THD_countmask( nvox , mask ) ;
647 INC_MEM_STATS( nmask * sizeof(byte), "mask array" );
648 PRINT_MEM_STATS( "mask" );
649
650 INFO_message("%d voxels in -mask dataset",nmask) ;
651 if( nmask < 2 ) ERROR_EXIT_CC("Only %d voxels in -mask, exiting...",nmask);
652
653 /* update running memory statistics to reflect loading the image */
654 DEC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );
655 PRINT_MEM_STATS( "mset unload" );
656
657 /* free all memory associated with the mask datast */
658 DSET_unload(mset) ;
659 DSET_delete(mset) ;
660 mset = NULL ;
661 }
662 /* if automasking is requested, handle that now */
663 else {
664 mask = THD_automask( xset ) ;
665 nmask = THD_countmask( nvox , mask ) ;
666 INC_MEM_STATS( nmask * sizeof(byte), "mask array" );
667 PRINT_MEM_STATS( "mask" );
668 INFO_message("No mask provided, created one using AFNI's automask procedure, %d voxels survive",nmask) ;
669 if( nmask < 2 ) ERROR_EXIT_CC("Only %d voxels in -automask!",nmask);
670 }
671
672 /*-- create vectim from input dataset --*/
673 INFO_message("vectim-izing input dataset") ;
674
675 /*-- CC added in mask to reduce the size of xvectim -- */
676 xvectim = THD_dset_to_vectim( xset , mask , 0 ) ;
677 if( xvectim == NULL ) ERROR_EXIT_CC("Can't create vectim?!") ;
678
679 /*-- CC update our memory stats to reflect vectim -- */
680 INC_MEM_STATS((xvectim->nvec*sizeof(int)) +
681 ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
682 sizeof(MRI_vectim), "vectim");
683 PRINT_MEM_STATS( "vectim" );
684
685 /* iterate through and remove all voxels that have zero
686 variance */
687 for (ii=0; ii<xvectim->nvec; ii++)
688 {
689 double sum = 0.0;
690 double sum_sq = 0.0;
691 int tt;
692
693 float* xsar = VECTIM_PTR(xvectim,ii);
694
695 for(tt=0; tt<xvectim->nvals; tt++)
696 {
697 sum += xsar[tt];
698 sum_sq += xsar[tt] * xsar[tt];
699 }
700
701 if((sum_sq-sum*sum/(double)nvals)/(double)(nvals-1) < eps)
702 {
703 mask[xvectim->ivec[ii]] = 0;
704 rem_count++;
705 }
706 }
707
708 /* update running memory statistics to reflect freeing the vectim */
709 DEC_MEM_STATS(((xvectim->nvec*sizeof(int)) +
710 ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
711 sizeof(MRI_vectim)), "vectim");
712
713 /* toss some trash */
714 VECTIM_destroy(xvectim) ;
715
716 /* make another xvectim with the updated mask */
717 xvectim = THD_dset_to_vectim( xset , mask , 0 ) ;
718 if( xvectim == NULL ) ERROR_EXIT_CC("Can't create vectim?!") ;
719
720 /*-- CC update our memory stats to reflect vectim -- */
721 INC_MEM_STATS((xvectim->nvec*sizeof(int)) +
722 ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
723 sizeof(MRI_vectim), "vectim");
724 PRINT_MEM_STATS( "vectim" );
725
726 INFO_message("!! %d voxels remain after removing those with no variance (%d)\n",
727 xvectim->nvec, rem_count);
728
729
730 /* unload DSET to reduce the amount of memory used */
731 DSET_unload(xset);
732 DEC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");
733
734 /** For the case of Pearson correlation, we make sure the **/
735 /** data time series have their mean removed (polort >= 0) **/
736 /** and are normalized, so that correlation = dot product, **/
737 /** and we can use function zm_THD_pearson_corr for speed. **/
738
739 switch( method ){
740 default:
741 case PEARSON: corfun = zm_THD_pearson_corr ; break ;
742 case ETA2: corfun = my_THD_eta_squared ; break ;
743 }
744
745 if( method == ETA2 && polort >= 0 )
746 WARNING_message("Polort for -eta2 should probably be -1...");
747
748
749 /*--- CC the vectim contains a mapping between voxel index and mask index,
750 tap into that here to avoid duplicating memory usage, also create a
751 mapping that goes the other way ---*/
752
753 if( mask != NULL )
754 {
755 /* tap into the xvectim mapping */
756 mask_ndx_to_vol_ndx = xvectim->ivec;
757
758 /* create a mapping that goes the opposite way */
759 if(( vol_ndx_to_mask_ndx = (long*)calloc( nvox, sizeof(long) )) == NULL)
760 {
761 ERROR_EXIT_CC(
762 "Could not allocate %d byte array for mask index mapping\n",
763 nmask*sizeof(long));
764 }
765
766 /* --- CC account for the size of the mask */
767 INC_MEM_STATS( nvox*sizeof(long), "ndx mask array" );
768
769 /* --- create the reverse mapping */
770 for ( ii = 0; ii < xvectim->nvec; ii++)
771 {
772 vol_ndx_to_mask_ndx[ mask_ndx_to_vol_ndx[ ii ] ] = ii;
773 }
774
775 /* --- CC free the mask */
776 DEC_MEM_STATS( nmask*sizeof(byte), "mask array" );
777 free(mask); mask=NULL;
778 PRINT_MEM_STATS( "mask unload" );
779 }
780
781 /* -- CC configure detrending --*/
782 if( polort < 0 && method == PEARSON )
783 {
784 polort = 0;
785 WARNING_message("Pearson correlation always uses polort >= 0");
786 }
787 if( polort >= 0 )
788 {
789 for( ii=0 ; ii < xvectim->nvec ; ii++ )
790 {
791 /* remove polynomial trend */
792 DETREND_polort(polort,nvals,VECTIM_PTR(xvectim,ii)) ;
793 }
794 }
795
796
797 /* -- this procedure does not change time series that
798 have zero variance -- */
799 if( method == PEARSON ) THD_vectim_normalize(xvectim) ; /* L2 norm = 1 */
800
801 /*-- create output dataset --*/
802 cset = EDIT_empty_copy( xset ) ;
803
804 /*-- configure the output dataset */
805 if( abuc )
806 {
807 EDIT_dset_items(
808 cset ,
809 ADN_prefix , prefix ,
810 ADN_nvals , nsubbriks , /* 2 subbricks */
811 ADN_ntt , 0 , /* no time axis */
812 ADN_type , HEAD_ANAT_TYPE ,
813 ADN_func_type , ANAT_BUCK_TYPE ,
814 ADN_datum_all , MRI_float ,
815 ADN_none ) ;
816 }
817 else
818 {
819 EDIT_dset_items( cset ,
820 ADN_prefix , prefix ,
821 ADN_nvals , nsubbriks , /* 2 subbricks */
822 ADN_ntt , nsubbriks , /* num times */
823 ADN_ttdel , 1.0 , /* fake TR */
824 ADN_nsl , 0 , /* no slice offsets */
825 ADN_type , HEAD_ANAT_TYPE ,
826 ADN_func_type , ANAT_EPI_TYPE ,
827 ADN_datum_all , MRI_float ,
828 ADN_none ) ;
829 }
830
831 /* add history information to the hearder */
832 tross_Make_History( "3dLFCD" , argc,argv , cset ) ;
833
834 INFO_message("creating output dataset in memory") ;
835
836 /* -- Configure the subbriks: Binary LFCD */
837 subbrik = 0;
838 EDIT_BRICK_TO_NOSTAT(cset,subbrik) ; /* stat params */
839 /* CC this sets the subbrik scaling factor, which we will probably want
840 to do again after we calculate the voxel values */
841 EDIT_BRICK_FACTOR(cset,subbrik,1.0) ; /* scale factor */
842
843 sprintf(str,"Binary LFCD") ;
844
845 EDIT_BRICK_LABEL(cset,subbrik,str) ;
846 EDIT_substitute_brick(cset,subbrik,MRI_float,NULL) ; /* make array */
847
848 /* copy measure data into the subbrik */
849 bodset = DSET_ARRAY(cset,subbrik);
850
851 /* -- Configure the subbriks: Weighted LFCD */
852 subbrik = 1;
853 EDIT_BRICK_TO_NOSTAT(cset,subbrik) ; /* stat params */
854 /* CC this sets the subbrik scaling factor, which we will probably want
855 to do again after we calculate the voxel values */
856 EDIT_BRICK_FACTOR(cset,subbrik,1.0) ; /* scale factor */
857
858 sprintf(str,"Weighted LFCD") ;
859
860 EDIT_BRICK_LABEL(cset,subbrik,str) ;
861 EDIT_substitute_brick(cset,subbrik,MRI_float,NULL) ; /* make array */
862
863 /* copy measure data into the subbrik */
864 wodset = DSET_ARRAY(cset,subbrik);
865
866 /* increment memory stats */
867 INC_MEM_STATS( (DSET_NVOX(cset)*DSET_NVALS(cset)*sizeof(float)),
868 "output dset");
869 PRINT_MEM_STATS( "outset" );
870
871 /*-- tell the user what we are about to do --*/
872 INFO_message( "Calculating LFCD with threshold = %f.\n", thresh);
873
874 /*---------- loop over mask voxels, correlate ----------*/
875 AFNI_OMP_START ;
876 #pragma omp parallel if( nmask > 999 )
877 {
878 int dx = 0;
879 int dy = 0;
880 int dz = 0;
881 int ix = 0;
882 int jy = 0;
883 int kz = 0;
884 int target_mask_ndx = 0;
885 int lout,ithr,nthr,vstep,vii ;
886 int neighbor_index;
887 float *xsar , *ysar ;
888 list_node* current_node = NULL ;
889 list_node* new_node = NULL ;
890 list_node* recycled_nodes = NULL ;
891 list_node* boundary_list = NULL ;
892 long* seen_voxels;
893 double car = 0.0 ;
894
895 /*-- get information about who we are --*/
896 #ifdef USE_OMP
897 ithr = omp_get_thread_num() ;
898 nthr = omp_get_num_threads() ;
899 if( ithr == 0 ) INFO_message("%d OpenMP threads started",nthr) ;
900 #else
901 ithr = 0 ; nthr = 1 ;
902 #endif
903
904 /*-- For the progress tracker, we want to print out 50 numbers,
905 figure out a number of loop iterations that will make this easy */
906 vstep = (int)( nmask / (nthr*50.0f) + 0.901f ) ; vii = 0 ;
907 if((MEM_STAT==0) && (ithr == 0 )) fprintf(stderr,"Looping:") ;
908
909 /* allocate a vector to track previously seen voxels */
910 if((seen_voxels = (long*)calloc(xvectim->nvec,sizeof(long)))==NULL)
911 {
912 ERROR_EXIT_CC( "Could not allocate memory for seen_voxels!\n");
913 }
914
915 #pragma omp for schedule(static, 1)
916 for( lout=0 ; lout < xvectim->nvec ; lout++ )
917 { /*----- outer voxel loop -----*/
918
919 if( ithr == 0 && vstep > 2 )
920 {
921 vii++ ;
922 if( vii%vstep == vstep/2 && MEM_STAT == 0 ) vstep_print();
923 }
924
925 /* get ref time series from this voxel */
926 xsar = VECTIM_PTR(xvectim,lout) ;
927
928 /* initialize the boundary with the target voxel */
929 if( recycled_nodes == NULL)
930 {
931 /* this looks like it is redundant, but I want the first if
932 statement to run in the critical section and the second
933 if statement to run after it */
934 #pragma omp critical(mem_alloc)
935 if(( new_node = (list_node*)malloc(sizeof(list_node)) ) != NULL )
936 {
937 INC_MEM_STATS((sizeof(list_node)), "list nodes");
938 }
939 if( new_node == NULL )
940 {
941 WARNING_message( "Could not allocate list node\n");
942 continue;
943 }
944 }
945 else
946 {
947 new_node = recycled_nodes;
948 recycled_nodes = recycled_nodes->next;
949 }
950
951 /* determine the full ndx for the seed target */
952 new_node->vox_vol_ndx = mask_ndx_to_vol_ndx[ lout ];
953
954 /* add source, dest, correlation to 1D file */
955 new_node->ix = DSET_index_to_ix(cset, new_node->vox_vol_ndx) ;
956 new_node->jy = DSET_index_to_jy(cset, new_node->vox_vol_ndx) ;
957 new_node->kz = DSET_index_to_kz(cset, new_node->vox_vol_ndx) ;
958
959 /* push the new node onto the boundary list,
960 if the boundary list is empty, there is a problem */
961 if (boundary_list != NULL)
962 {
963 WARNING_message("Boundary list not empty!\n");
964 }
965 new_node->next = boundary_list;
966 boundary_list = new_node;
967
968 /* reset the seen_voxels map */
969 memset(seen_voxels, 0, xvectim->nvec*sizeof(long));
970
971 /* iterate over the boundary until it is empty */
972 while( boundary_list != NULL )
973 {
974
975 /* pop a node off of the list */
976 current_node = boundary_list;
977 boundary_list = boundary_list->next;
978
979 /* iterate through a box around the current voxel */
980 for ( neighbor_index = 0; neighbor_index < num_neighbors; neighbor_index++ )
981 {
982 dx = neighborhood[3*neighbor_index+X];
983 dy = neighborhood[3*neighbor_index+Y];
984 dz = neighborhood[3*neighbor_index+Z];
985
986 ix = ( current_node->ix + (dx-1) );
987 /* make sure that we are in bounds */
988 if (( ix < 0 ) || ( ix > DSET_NX(xset) ))
989 {
990 continue;
991 }
992
993 jy = ( current_node->jy + (dy-1) );
994 /* make sure that we are in bounds */
995 if (( jy < 0 ) || ( jy > DSET_NY(xset) ))
996 {
997 continue;
998 }
999
1000 kz = ( current_node->kz + (dz-1) );
1001 /* make sure that we are in bounds */
1002 if (( kz < 0 ) || (kz > DSET_NZ(xset)))
1003 {
1004 continue;
1005 }
1006
1007 /* get the index of this voxel */
1008 target_mask_ndx =
1009 vol_ndx_to_mask_ndx[ DSET_ixyz_to_index(xset,ix,jy,kz) ];
1010
1011 /* if the voxel is in the mask, and hasn't been
1012 seen before, evaluate it for inclusion in the
1013 boundary */
1014 if(( target_mask_ndx != 0 ) && ( seen_voxels[ target_mask_ndx ] != 1 ))
1015 {
1016
1017 /* indicate that we have seen the voxel */
1018 seen_voxels[ target_mask_ndx ] = 1;
1019
1020 /* extract the time series */
1021 ysar = VECTIM_PTR(xvectim,target_mask_ndx) ;
1022
1023 /* calculate the correlation */
1024 car = (double)(corfun(nvals,xsar,ysar)) ;
1025
1026 /* if correlation is above threshold, add
1027 it to the LFCD stats and to the boundary */
1028 if ( car > thresh )
1029 {
1030
1031 if( recycled_nodes == NULL)
1032 {
1033 /* this looks like it is redundant, but I want the first if
1034 statement to run in the critical section and the second
1035 if statement to run after it */
1036 #pragma omp critical(mem_alloc)
1037 if(( new_node = (list_node*)malloc(sizeof(list_node)) ) != NULL )
1038 {
1039 INC_MEM_STATS((sizeof(list_node)), "list nodes");
1040 }
1041 if( new_node == NULL )
1042 {
1043 WARNING_message( "Could not allocate list node\n");
1044 continue;
1045 }
1046 }
1047 else
1048 {
1049 new_node = recycled_nodes;
1050 recycled_nodes = recycled_nodes->next;
1051 }
1052
1053 /* determine the full ndx for the seed target */
1054 new_node->vox_vol_ndx = mask_ndx_to_vol_ndx[ target_mask_ndx ];
1055
1056 /* add source, dest, correlation to 1D file */
1057 new_node->ix = ix ;
1058 new_node->jy = jy ;
1059 new_node->kz = kz ;
1060
1061 /* add the node to the boundary */
1062 new_node->next = boundary_list;
1063 boundary_list = new_node;
1064
1065 /* now increment the LFCD measures, this is
1066 done in a critical section */
1067 #pragma omp critical(dataupdate)
1068 {
1069 wodset[ mask_ndx_to_vol_ndx[ lout ]] += car;
1070 bodset[ mask_ndx_to_vol_ndx[ lout ]] += 1;
1071 }
1072 } /* if car > thresh */
1073 } /* if vox is in mask and hasn't been seen */
1074 } /* for neighbor_index */
1075
1076 /* we have finished processing this node, so recycle it */
1077 current_node->next = recycled_nodes;
1078 recycled_nodes = current_node;
1079
1080 } /* while( boundary_list != NULL ) */
1081 } /* for lout */
1082
1083 /* clean up the memory that is local to this thread */
1084 if (seen_voxels != NULL) free( seen_voxels );
1085
1086 /* clean out the recycled list */
1087 #pragma omp critical(mem_alloc)
1088 FREE_LIST( recycled_nodes )
1089 PRINT_MEM_STATS( "Free recycled nodes" );
1090
1091 } /* end OpenMP */
1092 AFNI_OMP_END ;
1093
1094 /* update the user so that they know what we are up to */
1095 INFO_message ("AFNI_OMP finished\n");
1096
1097 INFO_message("Done..\n") ;
1098
1099 if( vol_ndx_to_mask_ndx != NULL )
1100 {
1101 DEC_MEM_STATS(nvox*sizeof(long), "ndx mask array");
1102 free(vol_ndx_to_mask_ndx);
1103 vol_ndx_to_mask_ndx = NULL;
1104 PRINT_MEM_STATS( "Free ndx mask array" );
1105 }
1106
1107 /* update running memory statistics to reflect freeing the vectim */
1108 DEC_MEM_STATS(((xvectim->nvec*sizeof(int)) +
1109 ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
1110 sizeof(MRI_vectim)), "vectim");
1111
1112 /* toss some trash */
1113 VECTIM_destroy(xvectim) ;
1114
1115
1116 DSET_delete(xset) ;
1117
1118 PRINT_MEM_STATS( "Vectim and Input dset Unload" );
1119
1120 /* finito */
1121 INFO_message("Writing output dataset to disk [%s bytes]",
1122 commaized_integer_string(cset->dblk->total_bytes)) ;
1123
1124 /* write the dataset */
1125 DSET_write(cset) ;
1126 WROTE_DSET(cset) ;
1127
1128 /* increment our memory stats, since we are relying on the header for this
1129 information, we update the stats before actually freeing the memory */
1130 DEC_MEM_STATS( (DSET_NVOX(cset)*DSET_NVALS(cset)*sizeof(float)), "output dset");
1131
1132 /* free up the output dataset memory */
1133 DSET_unload(cset) ;
1134 DSET_delete(cset) ;
1135
1136 PRINT_MEM_STATS( "Unload Output Dset" );
1137
1138 /* force a print */
1139 MEM_STAT = 1;
1140 PRINT_MEM_STATS( "Fin" );
1141
1142 exit(0) ;
1143 }
1144