1 #ifdef USE_OMP
2 #include <omp.h>
3 #endif
4
5 #include "mrilib.h"
6 #include <sys/mman.h>
7 #include <sys/types.h>
8
9 #define SPEARMAN 1
10 #define QUADRANT 2
11 #define PEARSON 3
12 #define ETA2 4
13
14 #define ALLOW_MMAP
15
16 #ifndef ALLOW_MMAP
17 # define do_mmap 0
18 #else
19 static int do_mmap = 0 ;
20 #endif
21
22 /*----------------------------------------------------------------*/
23 /**** Include these here for potential optimization for OpenMP ****/
24 /*----------------------------------------------------------------*/
25 /*! Pearson correlation of x[] and y[] (x and y are NOT modified.
26 And we know ahead of time that the time series have 0 mean
27 and L2 norm 1.
28 *//*--------------------------------------------------------------*/
29
zm_THD_pearson_corr(int n,float * x,float * y)30 float zm_THD_pearson_corr( int n, float *x , float *y ) /* inputs are */
31 { /* zero mean */
32 register float xy ; register int ii ; /* and norm=1 */
33 if( n%2 == 0 ){
34 xy = 0.0f ;
35 for( ii=0 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
36 } else {
37 xy = x[0]*y[0] ;
38 for( ii=1 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
39 }
40 return xy ;
41 }
42
43 /*----------------------------------------------------------------*/
44 /* General correlation calculation. */
45
46 #if 0
47 float my_THD_pearson_corr( int n, float *x , float *y )
48 {
49 float xv,yv,xy , vv,ww , xm,ym ;
50 register int ii ;
51
52 xm = ym = 0.0f ;
53 for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
54 xm /= n ; ym /= n ;
55 xv = yv = xy = 0.0f ;
56 for( ii=0 ; ii < n ; ii++ ){
57 vv = x[ii]-xm ; ww = y[ii]-ym ; xv += vv*vv ; yv += ww*ww ; xy += vv*ww ;
58 }
59
60 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
61 return xy/sqrtf(xv*yv) ;
62 }
63 #endif
64
65 /*----------------------------------------------------------------*/
66 /*! eta^2 (Cohen, NeuroImage 2008) 25 Jun 2010 [rickr]
67 *
68 * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ]
69 * ------------------------------------
70 * SUM[ (a_i - M )^2 + (b_i - M )^2 ]
71 *
72 * where o a_i and b_i are the vector elements
73 * o m_i = (a_i + b_i)/2
74 * o M = mean across both vectors
75 -----------------------------------------------------------------*/
76
my_THD_eta_squared(int n,float * x,float * y)77 float my_THD_eta_squared( int n, float *x , float *y )
78 {
79 float num , denom , gm , lm, vv, ww;
80 register int ii ;
81
82 gm = 0.0f ;
83 for( ii=0 ; ii < n ; ii++ ){ gm += x[ii] + y[ii] ; }
84 gm /= (2.0f*n) ;
85
86 num = denom = 0.0f ;
87 for( ii=0 ; ii < n ; ii++ ){
88 lm = 0.5f * ( x[ii] + y[ii] ) ;
89 vv = (x[ii]-lm); ww = (y[ii]-lm); num += ( vv*vv + ww*ww );
90 vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww );
91 }
92
93 if( num < 0.0f || denom <= 0.0f || num >= denom ) return 0.0f ;
94 return (1.0f - num/denom) ;
95 }
96
97 /*----------------------------------------------------------------------------*/
98
vstep_print(void)99 static void vstep_print(void)
100 {
101 static int nn=0 ;
102 static char xx[10] = "0123456789" ;
103 fprintf(stderr , "%c" , xx[nn%10] ) ;
104 if( nn%10 == 9) fprintf(stderr,",") ;
105 nn++ ;
106 }
107
108 /*-----------------------------------------------------------------------------*/
109
110 #ifdef ALLOW_MMAP
111
112 static char *cbrik = NULL ;
113 static int64_t ncbrik = 0 ;
114 static char *cbrik_name = NULL ;
115
116 #undef MFLAG
117 #if defined(MAP_NOCACHE)
118 # define MFLAG (MAP_SHARED | MAP_NOCACHE)
119 #elif defined(MAP_NORESERVE)
120 # define MFLAG (MAP_SHARED | MAP_NORESERVE)
121 #else
122 # define MFLAG MAP_SHARED
123 #endif
124
125 #include <signal.h>
126
AC_sigfunc(int sig)127 void AC_sigfunc(int sig) /** signal handler for fatal errors **/
128 {
129 char *sname ;
130 static volatile int fff=0 ;
131 if( fff ) _exit(1) ; else fff = 1 ;
132 switch(sig){
133 default: sname = "unknown" ; break ;
134 case SIGINT: sname = "SIGINT" ; break ;
135 case SIGPIPE: sname = "SIGPIPE" ; break ;
136 case SIGSEGV: sname = "SIGSEGV" ; break ;
137 case SIGBUS: sname = "SIGBUS" ; break ;
138 case SIGTERM: sname = "SIGTERM" ; break ;
139 }
140 fprintf(stderr,"\n** 3dAutoTcorrelate: Fatal Signal %d (%s) received\n",sig,sname) ;
141 if( cbrik != NULL ){
142 fprintf(stderr,"** Un-mmap-ing %s file (but not deleting it)\n",cbrik_name) ;
143 munmap(cbrik,ncbrik) ; cbrik = NULL ;
144 }
145 exit(1) ;
146 }
147
148 #endif /* ALLOW_MMAP */
149
150 /*----------------------------------------------------------------------------*/
151
main(int argc,char * argv[])152 int main( int argc , char *argv[] )
153 {
154 THD_3dim_dataset *xset , *cset, *mset=NULL, *msetinner=NULL ;
155 int nopt=1 , method=PEARSON , do_autoclip=0 ;
156 int nvox , nvals , ii,kout , polort=1 , ix,jy,kz ;
157 char *prefix = "ATcorr" ;
158 byte *mask=NULL, *maskinner=NULL;
159 int nmask , nmaskinner , abuc=1 ;
160 int all_source=0; /* output all source voxels 25 Jun 2010 [rickr] */
161 char str[32] , *cpt ;
162 int *imap ; MRI_vectim *xvectim ;
163 float (*corfun)(int,float *,float*) = NULL ;
164
165 FILE *fout1D=NULL;
166
167 /*----*/
168
169 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
170
171 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
172 printf(
173 "Usage: 3dAutoTcorrelate [options] dset\n"
174 "Computes the correlation coefficient between the time series of each\n"
175 "pair of voxels in the input dataset, and stores the output into a\n"
176 "new anatomical bucket dataset [scaled to shorts to save memory space].\n"
177 "\n"
178 "*** Also see program 3dTcorrMap ***\n"
179 "\n"
180 "Options:\n"
181 " -pearson = Correlation is the normal Pearson (product moment)\n"
182 " correlation coefficient [default].\n"
183 " -eta2 = Output is eta^2 measure from Cohen et al., NeuroImage, 2008:\n"
184 " http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2705206/\n"
185 " http://dx.doi.org/10.1016/j.neuroimage.2008.01.066\n"
186 " ** '-eta2' is intended to be used to measure the similarity\n"
187 " between 2 correlation maps; therefore, this option is\n"
188 " to be used in a second stage analysis, where the input\n"
189 " dataset is the output of running 3dAutoTcorrelate with\n"
190 " the '-pearson' option -- the voxel 'time series' from\n"
191 " that first stage run is the correlation map of that\n"
192 " voxel with all other voxels.\n"
193 " ** '-polort -1' is recommended with this option!\n"
194 " ** Odds are you do not want use this option if the dataset\n"
195 " on which eta^2 is to be computed was generated with\n"
196 " options -mask_only_targets or -mask_source.\n"
197 " In this program, the eta^2 is computed between pseudo-\n"
198 " timeseries (the 4th dimension of the dataset).\n"
199 " If you want to compute eta^2 between sub-bricks then use\n"
200 " 3ddot -eta2 instead.\n"
201 #if 0
202 " -spearman = Correlation is the Spearman (rank) correlation\n"
203 " coefficient.\n"
204 " -quadrant = Correlation is the quadrant correlation coefficient.\n"
205 #else
206 " -spearman AND -quadrant are disabled at this time :-(\n"
207 #endif
208 "\n"
209 " -polort m = Remove polynomial trend of order 'm', for m=-1..3.\n"
210 " [default is m=1; removal is by least squares].\n"
211 " Using m=-1 means no detrending; this is only useful\n"
212 " for data/information that has been pre-processed.\n"
213 "\n"
214 " -autoclip = Clip off low-intensity regions in the dataset,\n"
215 " -automask = so that the correlation is only computed between\n"
216 " high-intensity (presumably brain) voxels. The\n"
217 " mask is determined the same way that 3dAutomask works.\n"
218 "\n"
219 " -mask mmm = Mask of both 'source' and 'target' voxels.\n"
220 " ** Restricts computations to those in the mask. Output\n"
221 " volumes are restricted to masked voxels. Also, only\n"
222 " masked voxels will have non-zero output.\n"
223 " ** A dataset with 1000 voxels would lead to output of\n"
224 " 1000 sub-bricks. With a '-mask' of 50 voxels, the\n"
225 " output dataset have 50 sub-bricks, where the 950\n"
226 " unmasked voxels would be all zero in all 50 sub-bricks\n"
227 " (unless option '-mask_only_targets' is also used).\n"
228 " ** The mask is encoded in the output dataset header in the\n"
229 " attribute named 'AFNI_AUTOTCORR_MASK' (cf. 3dMaskToASCII).\n"
230 "\n"
231 " -mask_only_targets = Provide output for all voxels.\n"
232 " ** Used with '-mask': every voxel is correlated with each\n"
233 " of the mask voxels. In the example above, there would\n"
234 " be 50 output sub-bricks; the n-th output sub-brick\n"
235 " would contain the correlations of the n-th voxel in\n"
236 " the mask with ALL 1000 voxels in the dataset (rather\n"
237 " than with just the 50 voxels in the mask).\n"
238 "\n"
239 " -mask_source sss = Provide ouput for voxels only in mask sss\n"
240 " ** For each seed in mask mm, compute correlations only with \n"
241 " non-zero voxels in sss. If you have 250 non-zero voxels \n"
242 " in sss, then the output will still have 50 sub-bricks, but\n"
243 " each n-th sub-brick will have non-zero values at the 250\n"
244 " non-zero voxels in sss\n"
245 " Do not use this option along with -mask_only_targets\n"
246 "\n"
247 " -prefix p = Save output into dataset with prefix 'p'\n"
248 " [default prefix is 'ATcorr'].\n"
249 " -out1D FILE.1D = Save output in a text file formatted thusly:\n"
250 " Row 1 contains the 1D indices of non zero voxels in the \n"
251 " mask from option -mask.\n"
252 " Column 1 contains the 1D indices of non zero voxels in the\n"
253 " mask from option -mask_source\n"
254 " The rest of the matrix contains the correlation/eta2 \n"
255 " values. Each column k corresponds to sub-brick k in \n"
256 " the output volume p.\n"
257 " To see 1D indices in AFNI, right click on the top left\n"
258 " corner of the AFNI controller - where coordinates are\n"
259 " shown - and chose voxel indices.\n"
260 " A 1D index (ijk) is computed from the 3D (i,j,k) indices:\n"
261 " ijk = i + j*Ni + k*Ni*Nj , with Ni and Nj being the\n"
262 " number of voxels in the slice orientation and given by:\n"
263 " 3dinfo -ni -nj YOUR_VOLUME_HERE\n"
264 " This option can only be used in conjunction with \n"
265 " options -mask and -mask_source. Otherwise it makes little\n"
266 " sense to write a potentially enormous text file.\n"
267 "\n"
268 " -time = Mark output as a 3D+time dataset instead of an anat bucket.\n"
269 #ifdef ALLOW_MMAP
270 "\n"
271 " -mmap = Write .BRIK results to disk directly using Unix mmap().\n"
272 " This trick can speed the program up when the amount\n"
273 " of memory required to hold the output is very large.\n"
274 " ** In many case, the amount of time needed to write\n"
275 " the results to disk is longer than the CPU time.\n"
276 " This option can shorten the disk write time.\n"
277 " ** If the program crashes, you'll have to manually\n"
278 " remove the .BRIK file, which will have been created\n"
279 " before the loop over voxels and written into during\n"
280 " that loop, rather than being written all at once\n"
281 " at the end of the analysis, as is usually the case.\n"
282 " ** If the amount of memory needed is bigger than the\n"
283 " RAM on your system, this program will be very slow\n"
284 " with or without '-mmap'.\n"
285 " ** This option won't work with NIfTI-1 (.nii) output!\n"
286 #else
287 "\n"
288 " -mmap is disabled at this time :-(\n"
289 #endif
290 "\n"
291 "Example: correlate every voxel in mask_in+tlrc with only those voxels in\n"
292 " mask_out+tlrc (the rest of each volume is zero, for speed).\n"
293 " Assume detrending was already done along with other pre-processing.\n"
294 " The output will have one volume per masked voxel in mask_in+tlrc.\n"
295 " Volumes will be labeled by the ijk index triples of mask_in+tlrc.\n"
296 "\n"
297 " 3dAutoTcorrelate -mask_source mask_out+tlrc -mask mask_in+tlrc \\\n"
298 " -polort -1 -prefix test_corr clean_epi+tlrc\n"
299 "\n"
300 "Notes:\n"
301 " * The output dataset is anatomical bucket type of shorts\n"
302 " (unless '-time' is used).\n"
303 " * Values are scaled so that a correlation (or eta-squared)\n"
304 " of 1 corresponds to a value of 10000.\n"
305 " * The output file might be gigantic and you might run out\n"
306 " of memory running this program. Use at your own risk!\n"
307 " ++ If you get an error message like\n"
308 " *** malloc error for dataset sub-brick\n"
309 " this means that the program ran out of memory when making\n"
310 " the output dataset.\n"
311 #ifdef ALLOW_MMAP
312 " ++ If this happens, you can try to use the '-mmap' option,\n"
313 " and if you are lucky, the program may actually run.\n"
314 #endif
315 " * The program prints out an estimate of its memory usage\n"
316 " when it starts. It also prints out a progress 'meter'\n"
317 " to keep you pacified.\n"
318 " * This is a quick hack for Peter Bandettini. Now pay up.\n"
319 " * OpenMP-ized for Hang Joon Jo. Where's my baem-sul?\n"
320 "\n"
321 "-- RWCox - 31 Jan 2002 and 16 Jul 2010\n"
322 ) ;
323 PRINT_AFNI_OMP_USAGE("3dAutoTcorrelate",NULL) ;
324 PRINT_COMPILE_DATE ; exit(0) ;
325 }
326
327 mainENTRY("3dAutoTcorrelate main"); machdep(); PRINT_VERSION("3dAutoTcorrelate");
328 AFNI_logger("3dAutoTcorrelate",argc,argv);
329
330 /*-- option processing --*/
331
332 while( nopt < argc && argv[nopt][0] == '-' ){
333
334 #ifdef ALLOW_MMAP
335 if( strcmp(argv[nopt],"-mmap") == 0 ){ /* Jul 2010 */
336 do_mmap = 1 ; nopt++ ; continue ;
337 }
338 #endif
339
340 if( strcmp(argv[nopt],"-time") == 0 ){
341 abuc = 0 ; nopt++ ; continue ;
342 }
343
344 if( strcmp(argv[nopt],"-autoclip") == 0 ||
345 strcmp(argv[nopt],"-automask") == 0 ){
346
347 do_autoclip = 1 ; nopt++ ; continue ;
348 }
349
350 if( strcmp(argv[nopt],"-mask") == 0 ){
351 mset = THD_open_dataset(argv[++nopt]);
352 CHECK_OPEN_ERROR(mset,argv[nopt]);
353 nopt++ ; continue ;
354 }
355
356 if( strcmp(argv[nopt],"-mask_source") == 0 ){
357 msetinner = THD_open_dataset(argv[++nopt]);
358 CHECK_OPEN_ERROR(msetinner,argv[nopt]);
359 nopt++ ; continue ;
360 }
361
362 if( strcmp(argv[nopt],"-out1D") == 0 ){
363 if (!(fout1D = fopen(argv[++nopt], "w"))) {
364 ERROR_message("Failed to open %s for writing", argv[nopt]);
365 exit(1);
366 }
367 nopt++ ; continue ;
368 }
369
370 if( strcmp(argv[nopt],"-mask_only_targets") == 0 ){
371 all_source = 1 ; nopt++ ; continue ;
372 }
373
374 if( strcmp(argv[nopt],"-pearson") == 0 ){
375 method = PEARSON ; nopt++ ; continue ;
376 }
377
378 #if 0
379 if( strcmp(argv[nopt],"-spearman") == 0 ){
380 method = SPEARMAN ; nopt++ ; continue ;
381 }
382
383 if( strcmp(argv[nopt],"-quadrant") == 0 ){
384 method = QUADRANT ; nopt++ ; continue ;
385 }
386 #endif
387
388 if( strcmp(argv[nopt],"-eta2") == 0 ){
389 method = ETA2 ; nopt++ ; continue ;
390 }
391
392 if( strcmp(argv[nopt],"-prefix") == 0 ){
393 prefix = strdup(argv[++nopt]) ;
394 if( !THD_filename_ok(prefix) ){
395 ERROR_exit("Illegal value after -prefix!") ;
396 }
397 nopt++ ; continue ;
398 }
399
400 if( strcmp(argv[nopt],"-polort") == 0 ){
401 int val = (int)strtod(argv[++nopt],&cpt) ;
402 if( *cpt != '\0' || val < -1 || val > 3 ){
403 ERROR_exit("Illegal value after -polort!") ;
404 }
405 polort = val ; nopt++ ; continue ;
406 }
407
408 ERROR_exit("Illegal option: %s",argv[nopt]) ;
409 }
410
411 #ifdef ALLOW_MMAP
412 cpt = strstr(prefix,".nii") ; /* Jul 2010 */
413 if( do_mmap && cpt != NULL ){
414 *cpt = '\0' ; if( *prefix == '\0' ) prefix = "ATcorr" ;
415 WARNING_message("-mmap ==> can't use NIfTI-1 ==> output prefix = %s",prefix) ;
416 }
417 #endif
418
419 /*-- open dataset, check for legality --*/
420
421 if( nopt >= argc ) ERROR_exit("Need a dataset on command line!?") ;
422
423 if (fout1D && (!mset || !msetinner)) {
424 fclose(fout1D);
425 ERROR_message("Must use -mask and -mask_source with -out1D");
426 exit(1);
427 }
428
429 xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);
430 if( DSET_NVALS(xset) < 3 )
431 ERROR_exit("Input dataset %s does not have 3 or more sub-bricks!",argv[nopt]) ;
432 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
433
434 /*-- compute mask array, if desired --*/
435
436 nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
437
438 if( mset ){
439 if( DSET_NVOX(mset) != nvox )
440 ERROR_exit("Input and mask dataset differ in number of voxels!") ;
441 mask = THD_makemask(mset, 0, 1.0, 0.0) ;
442 nmask = THD_countmask( nvox , mask ) ;
443 INFO_message("%d voxels in -mask dataset",nmask) ;
444 if( nmask < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmask);
445 DSET_unload(mset) ;
446 } else if( do_autoclip ){
447 mask = THD_automask( xset ) ;
448 nmask = THD_countmask( nvox , mask ) ;
449 INFO_message("%d voxels survive -autoclip",nmask) ;
450 if( nmask < 2 ) ERROR_exit("Only %d voxels in -automask!",nmask);
451 } else {
452 nmask = nvox ;
453 INFO_message("computing for all %d voxels",nmask) ;
454 }
455 if (msetinner) {
456 if( DSET_NVOX(msetinner) != nvox )
457 ERROR_exit("Input and mask_source dataset differ in number of voxels!");
458 maskinner = THD_makemask(msetinner, 0, 1.0, 0.0) ;
459 nmaskinner = THD_countmask( nvox , maskinner ) ;
460 INFO_message("%d voxels in -mask dataset",nmaskinner) ;
461 if( nmaskinner < 2 )
462 ERROR_exit("Only %d voxels in -mask, exiting...",nmaskinner);
463 DSET_unload(msetinner) ;
464 DSET_delete(msetinner); msetinner = NULL;
465 } else {
466 if (mask) {
467 maskinner = mask;
468 nmaskinner = -nvox; /* a flag to be sure we remember this deed */
469 }
470 }
471
472 if( method == ETA2 && polort >= 0 )
473 WARNING_message("Polort for -eta2 should probably be -1...");
474
475 /*-- create vectim from input dataset --*/
476
477 /** For the case of Pearson correlation, we make sure the **/
478 /** data time series have their mean removed (polort >= 0) **/
479 /** and are normalized, so that correlation = dot product, **/
480 /** and we can use function zm_THD_pearson_corr for speed. **/
481
482 switch( method ){
483 default:
484 case PEARSON: corfun = zm_THD_pearson_corr ; break ;
485 case ETA2: corfun = my_THD_eta_squared ; break ;
486 }
487
488 INFO_message("vectim-izing input dataset") ;
489 xvectim = THD_dset_to_vectim( xset , NULL , 0 ) ;
490 if( xvectim == NULL ) ERROR_exit("Can't create vectim?!") ;
491 DSET_unload(xset) ;
492 if( polort < 0 && method == PEARSON ){
493 polort = 0; WARNING_message("Pearson correlation always uses polort >= 0");
494 }
495 if( polort >= 0 ){
496 for( ii=0 ; ii < nvox ; ii++ ){ /* remove polynomial trend */
497 DETREND_polort(polort,nvals,VECTIM_PTR(xvectim,ii)) ;
498 }
499 }
500 if( method == PEARSON ) THD_vectim_normalize(xvectim) ; /* L2 norm = 1 */
501
502 /*-- create output dataset --*/
503
504 cset = EDIT_empty_copy( xset ) ;
505
506 if( abuc ){
507 EDIT_dset_items( cset ,
508 ADN_prefix , prefix ,
509 ADN_nvals , nmask ,
510 ADN_ntt , 0 , /* no time axis */
511 ADN_type , HEAD_ANAT_TYPE ,
512 ADN_func_type , ANAT_BUCK_TYPE ,
513 ADN_datum_all , MRI_short ,
514 ADN_none ) ;
515 } else {
516 EDIT_dset_items( cset ,
517 ADN_prefix , prefix ,
518 ADN_nvals , nmask ,
519 ADN_ntt , nmask , /* num times */
520 ADN_ttdel , 1.0 , /* fake TR */
521 ADN_nsl , 0 , /* no slice offsets */
522 ADN_type , HEAD_ANAT_TYPE ,
523 ADN_func_type , ANAT_EPI_TYPE ,
524 ADN_datum_all , MRI_short ,
525 ADN_none ) ;
526 }
527
528 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(cset)) )
529 ERROR_exit("Output dataset %s already exists!",DSET_HEADNAME(cset)) ;
530
531 { double nb = (double)(xset->dblk->total_bytes)
532 +(double)(cset->dblk->total_bytes) ;
533 nb /= 1000000 ;
534 INFO_message(
535 "Memory required = %.1f Mbytes for %d output sub-bricks",nb,nmask);
536 if( nb > 2000.0 && (sizeof(void *) < 8 || sizeof(size_t) < 8) )
537 ERROR_exit("Can't run on a 32-bit system!") ;
538 #ifdef ALLOW_MMAP
539 if( nb > 1000.0 && !do_mmap )
540 INFO_message("If you run out of memory, try using the -mmap option.") ;
541 #endif
542 }
543
544 tross_Make_History( "3dAutoTcorrelate" , argc,argv , cset ) ;
545
546 /*--- make output dataset first ---*/
547
548 imap = (int *)calloc(sizeof(int),nmask) ;
549
550 if( !do_mmap )
551 ININFO_message("creating output dataset in memory") ;
552
553 for( kout=ii=0 ; ii < nvox ; ii++ ){
554 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip it */
555
556 EDIT_BRICK_TO_NOSTAT(cset,kout) ; /* stat params */
557 EDIT_BRICK_FACTOR(cset,kout,0.0001) ; /* scale factor */
558
559 ix = DSET_index_to_ix(cset,ii) ; /* brick label */
560 jy = DSET_index_to_jy(cset,ii) ;
561 kz = DSET_index_to_kz(cset,ii) ;
562 sprintf(str,"%03d_%03d_%03d",ix,jy,kz) ;
563 EDIT_BRICK_LABEL(cset,kout,str) ;
564 if( !do_mmap )
565 EDIT_substitute_brick(cset,kout,MRI_short,NULL) ; /* make array */
566
567 imap[kout] = ii ; kout++ ;
568 }
569
570 #ifdef ALLOW_MMAP
571 if( do_mmap ){
572 FILE *fp ; int64_t qq,dq ; char buf[2]={0,0} ;
573
574 signal(SIGINT ,AC_sigfunc) ; /* setup signal handler */
575 signal(SIGBUS ,AC_sigfunc) ; /* for fatal errors */
576 signal(SIGSEGV,AC_sigfunc) ;
577 signal(SIGTERM,AC_sigfunc) ;
578
579 cbrik_name = strdup(DSET_BRIKNAME(cset)) ; /* save name for later use */
580 ncbrik = cset->dblk->total_bytes ;
581 ININFO_message("creating output file %s via mmap() [%s bytes -- about %s]",
582 cbrik_name , commaized_integer_string(ncbrik) ,
583 approximate_number_string((double)ncbrik) ) ;
584 fp = fopen(cbrik_name,"w+") ;
585 if( fp == NULL ) ERROR_exit("Can't create output BRIK file :-(") ;
586 fseek( fp , ncbrik-1 , SEEK_SET ) ;
587 fwrite( &buf , 1 , 2 , fp ) ; fflush(fp) ;
588 cbrik = mmap( 0 , (size_t)ncbrik ,
589 PROT_READ | PROT_WRITE , MFLAG , fileno(fp) , 0 ) ;
590 fclose(fp) ;
591 if( cbrik == (char *)(-1) ){
592 remove(cbrik_name) ; ERROR_exit("Can't mmap output BRIK file :-(") ;
593 }
594 ININFO_message("Testing output BRIK") ;
595 dq = ncbrik / 65536 ; if( dq < 1024 ) dq = 1024 ;
596 for( qq=0 ; qq < ncbrik ; qq+=dq ) cbrik[qq] = 0 ;
597 cbrik[ncbrik-1] = 0 ;
598 ININFO_message("... done with test") ;
599 }
600 #endif
601
602 /*---------- loop over mask voxels, correlate ----------*/
603
604 AFNI_OMP_START ;
605 #pragma omp parallel if( nmask > 999 )
606 {
607 int ii,jj,kout , ithr,nthr , vstep,vii ;
608 float *xsar , *ysar ; short *car=NULL ;
609
610 #ifdef USE_OMP
611 ithr = omp_get_thread_num() ;
612 nthr = omp_get_num_threads() ;
613 if( ithr == 0 ) INFO_message("%d OpenMP threads started",nthr) ;
614 #else
615 ithr = 0 ; nthr = 1 ;
616 #endif
617
618 #ifdef ALLOW_MMAP
619 if( do_mmap ) car = (short *)malloc(sizeof(short)*nvox) ;
620 #endif
621
622 vstep = (int)( nmask / (nthr*50.0f) + 0.901f ) ; vii = 0 ;
623 if( ithr == 0 ) fprintf(stderr,"Looping:") ;
624
625 #pragma omp for
626 for( kout=0 ; kout < nmask ; kout++ ){ /*----- outer voxel loop -----*/
627
628 ii = imap[kout] ; /* ii= source voxel (we know that ii is in the mask) */
629
630 if( ithr == 0 && vstep > 2 ) /* allow small dsets 16 Jun 2011 [rickr] */
631 { vii++ ; if( vii%vstep == vstep/2 ) vstep_print(); }
632
633 /* get ref time series from this voxel */
634
635 xsar = VECTIM_PTR(xvectim,ii) ;
636 if( !do_mmap ) car = DSET_ARRAY(cset,kout) ;
637
638 for( jj=0 ; jj < nvox ; jj++ ){ /*----- inner loop over voxels -----*/
639
640 car[jj] = 0 ; /* default value */
641
642 /* skip unmasked voxels, unless want correlations from all source voxels */
643
644 if( !all_source && maskinner != NULL && maskinner[jj] == 0 ) continue ;
645
646 #if 0
647 if( jj == kout ){ car[jj] = 10000 ; continue ; } /* correlation = 1.0 */
648 #endif
649
650 ysar = VECTIM_PTR(xvectim,jj) ;
651
652 car[jj] = (short)(10000.49f*corfun(nvals,xsar,ysar)) ;
653
654 } /* end of inner loop over voxels */
655
656 #ifdef ALLOW_MMAP
657 if( do_mmap ){ /* copy results to disk mmap now */
658 short *cout = ((short *)(cbrik)) + (int64_t)(nvox)*(int64_t)(kout) ;
659 AAmemcpy( cout , car , sizeof(short)*nvox ) ;
660 }
661 #endif
662
663 } /* end of outer loop over ref voxels */
664
665 if( do_mmap ) free(car) ; /* workspace for each thread */
666 if( ithr == 0 ) fprintf(stderr,".") ;
667
668 } /* end OpenMP */
669 AFNI_OMP_END ;
670
671 /*---------- Finish up ---------*/
672
673 /* write mask info (if any) to output dataset header */
674
675 if( mask != NULL ){
676 char *maskstring = mask_to_b64string(nvox,mask) ;
677 THD_set_string_atr( cset->dblk , "AFNI_AUTOTCORR_MASK" , maskstring ) ;
678 free(maskstring) ; free(mask) ;
679 if( ISVALID_DSET(mset) ){
680 THD_set_string_atr( cset->dblk , "AFNI_AUTOCORR_MASK_IDCODE" ,
681 DSET_IDCODE_STR(mset) ) ;
682 THD_set_string_atr( cset->dblk , "AFNI_AUTOCORR_MASK_NAME" ,
683 DSET_HEADNAME(mset) ) ;
684 DSET_delete(mset) ;
685 }
686 }
687
688
689 /* toss some trash */
690
691 VECTIM_destroy(xvectim) ;
692
693 /* if did mmap(), finish that off as well */
694
695 if( !do_mmap ){
696 fprintf(stderr,"Done..\n") ;
697 } else {
698 #ifdef ALLOW_MMAP
699 char *ptr = cbrik ;
700 fprintf(stderr,"Done.[un-mmap-ing") ;
701 for( ii=0 ; ii < nmask ; ii++ ){ /* map sub-bricks */
702 mri_fix_data_pointer( ptr , DSET_BRICK(cset,ii) ) ;
703 ptr += DSET_BRICK_BYTES(cset,ii) ;
704 }
705 THD_load_statistics(cset) ; /* so can do stats */
706 fprintf(stderr,".") ;
707 munmap(cbrik,ncbrik) ; cbrik = NULL ;
708 fprintf(stderr,"].\n") ;
709 #endif
710 }
711
712 if (fout1D) { /* write results to ASCII file, hopefully when masked */
713 if (fout1D && (!mask || !maskinner)) {
714 ERROR_message("Option -1Dout restricted to commands using"
715 "-mask and -mask_source options.");
716 } else {
717 float *far = (float *)calloc(nmask, sizeof(float));
718 int jj;
719 fprintf(fout1D,"#Text output of:\n#");
720 for (ii=0; ii<argc; ++ii) fprintf(fout1D,"%s ", argv[ii]);
721 fprintf(fout1D,"\n"
722 "#First row contains 1D indices of voxels from -mask\n"
723 "#First column contains 1D indices of voxels from -mask_source\n");
724 /* write the matrix */
725 fprintf(fout1D," ");
726 for( ii=0 ; ii < nmask ; ii++ ) {
727 fprintf(fout1D,"%08d ", imap[ii]);
728 }
729 fprintf(fout1D," #Mask Voxel Indices (from -mask)");
730 for( jj=0 ; jj < nvox ; jj++ ){
731 if (maskinner[jj]) {
732 fprintf(fout1D,"\n%08d ", jj);
733 THD_extract_float_array( jj, cset, far );
734 for (ii=0; ii < nmask; ++ii)
735 fprintf(fout1D,"%f ", far[ii]);
736 }
737 }
738 free(far); far=NULL;
739 fclose(fout1D); fout1D = NULL;
740 }
741 }
742
743 /* free inner mask if not a copy of mask */
744 if ( maskinner && maskinner != mask) {
745 free(maskinner);
746 }
747 nmaskinner = 0;
748 maskinner=NULL;
749 if (imap) free(imap) ; imap = NULL;
750
751 /* finito */
752
753 if( !do_mmap ){
754 INFO_message("Writing output dataset to disk [%s bytes]",
755 commaized_integer_string(cset->dblk->total_bytes)) ;
756 THD_set_write_compression(COMPRESS_NONE) ; AFNI_setenv("AFNI_AUTOGZIP=NO") ;
757 DSET_write(cset) ;
758 } else {
759 INFO_message("Writing output dataset header") ;
760 DSET_write_header(cset) ;
761 }
762 WROTE_DSET(cset) ;
763
764 exit(0) ;
765 }
766