1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 #ifdef USE_OMP
8 #include "thd_Tcorr1D.c"
9 #endif
10
main(int argc,char * argv[])11 int main( int argc , char *argv[] )
12 {
13 THD_3dim_dataset *xset=NULL , *cset ;
14 int nopt=1, datum=MRI_float, nvals, ii;
15 MRI_IMAGE *ysim=NULL ;
16 char *prefix = "Tcorr1D", *smethod="pearson";
17 char *xnam=NULL , *ynam=NULL ;
18 byte *mask=NULL ; int mask_nx,mask_ny,mask_nz , nmask=0 ;
19 int do_atanh = 0 ; /* 12 Jan 2018 */
20
21 /*----*/
22
23 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
24
25 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
26 printf("\n"
27 "Usage: 3dTcorr1D [options] xset y1D ~1~\n"
28 "\n"
29 "Computes the correlation coefficient between each voxel time series\n"
30 "in the input 3D+time dataset 'xset' and each column in the 1D time\n"
31 "series file 'y1D', and stores the output values in a new dataset.\n"
32 "\n"
33 "--------\n"
34 "OPTIONS: ~1~\n"
35 "--------\n"
36 " -pearson = Correlation is the normal Pearson (product moment)\n"
37 " correlation coefficient [this is the default method].\n"
38 " -spearman = Correlation is the Spearman (rank) correlation\n"
39 " coefficient.\n"
40 " -quadrant = Correlation is the quadrant correlation coefficient.\n"
41 " -ktaub = Correlation is Kendall's tau_b coefficient.\n"
42 " ++ For 'continuous' or finely-discretized data, tau_b and\n"
43 " rank correlation are nearly equivalent (but not equal).\n"
44 " -dot = Doesn't actually compute a correlation coefficient; just\n"
45 " calculates the dot product between the y1D vector(s)\n"
46 " and the dataset time series.\n"
47 "\n"
48 " -Fisher = Apply the 'Fisher' (inverse hyperbolic tangent = arctanh)\n"
49 " transformation to the results.\n"
50 " ++ It does NOT make sense to use this with '-ktaub', but if\n"
51 " you want to do it, the program will not stop you.\n"
52 " ++ Cannot be used with '-dot'!\n"
53 "\n"
54 " -prefix p = Save output into dataset with prefix 'p'\n"
55 " [default prefix is 'Tcorr1D'].\n"
56 "\n"
57 " -mask mmm = Only process voxels from 'xset' that are nonzero\n"
58 " in the 3D mask dataset 'mmm'.\n"
59 " ++ Other voxels in the output will be set to zero.\n"
60 "\n"
61 " -float = Save results in float format [the default format].\n"
62 " -short = Save results in scaled short format [to save disk space].\n"
63 " ++ Cannot be used with '-dot'!\n"
64 "\n"
65 "------\n"
66 "NOTES: ~1~\n"
67 "------\n"
68 "\n"
69 "* The output dataset is functional bucket type, with one sub-brick\n"
70 " per column of the input y1D file.\n"
71 "\n"
72 "* No detrending, blurring, or other pre-processing options are available;\n"
73 " if you want these things, see 3dDetrend or 3dTproject or 3dcalc.\n"
74 " [In other words, this program presumes you know what you are doing!]\n"
75 "\n"
76 "* Also see 3dTcorrelate to do voxel-by-voxel correlation of TWO\n"
77 " 3D+time datasets' time series, with similar options.\n"
78 "\n"
79 "* You can extract the time series from a single voxel with given\n"
80 " spatial indexes using 3dmaskave, and then run it with 3dTcorr1D:\n"
81 " 3dmaskave -quiet -ibox 40 30 20 epi_r1+orig > r1_40_30_20.1D\n"
82 " 3dTcorr1D -pearson -Fisher -prefix c_40_30_20 epi_r1+orig r1_40_30_20.1D\n"
83 "\n"
84 "* http://en.wikipedia.org/wiki/Correlation\n"
85 "* http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient\n"
86 "* http://en.wikipedia.org/wiki/Spearman%%27s_rank_correlation_coefficient\n"
87 "* http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient\n"
88 "\n"
89 "-- RWCox - Apr 2010\n"
90 " - Jun 2010: Multiple y1D columns; OpenMP; -short; -mask.\n"
91 ) ;
92 PRINT_AFNI_OMP_USAGE("3dTcorr1D",NULL) ;
93 PRINT_COMPILE_DATE ; exit(0) ;
94 }
95
96 mainENTRY("3dTcorr1D main"); machdep(); AFNI_logger("3dTcorr1D",argc,argv);
97 PRINT_VERSION("3dTcorr1D") ; THD_check_AFNI_version("3dTcorr1D") ;
98
99 /*-- option processing --*/
100
101 while( nopt < argc && argv[nopt][0] == '-' ){
102
103 if( strcasecmp(argv[nopt],"-mask") == 0 ){ /* 28 Jun 2010 */
104 THD_3dim_dataset *mset ;
105 if( ++nopt >= argc ) ERROR_exit("Need argument after '-mask'") ;
106 if( mask != NULL ) ERROR_exit("Can't have two mask inputs") ;
107 mset = THD_open_dataset( argv[nopt] ) ;
108 CHECK_OPEN_ERROR(mset,argv[nopt]) ;
109 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
110 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
111 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
112 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[nopt]) ;
113 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
114 INFO_message("Number of voxels in mask = %d",nmask) ;
115 if( nmask < 2 ) ERROR_exit("Mask is too small to process") ;
116 nopt++ ; continue ;
117 }
118
119 if( strcasecmp(argv[nopt],"-float") == 0 ){ /* 27 Jun 2010 */
120 datum = MRI_float ; nopt++ ; continue ;
121 }
122 if( strcasecmp(argv[nopt],"-short") == 0 ){
123 datum = MRI_short ; nopt++ ; continue ;
124 }
125
126 if( strcasecmp(argv[nopt],"-pearson") == 0 ){
127 smethod = "pearson" ; nopt++ ; continue ;
128 }
129
130 if( strcasecmp(argv[nopt],"-dot") == 0 ){
131 smethod = "dot" ; nopt++ ; continue ;
132 }
133
134 if( strcasecmp(argv[nopt],"-spearman") == 0 || strcasecmp(argv[nopt],"-rank") == 0 ){
135 smethod = "spearman" ; nopt++ ; continue ;
136 }
137
138 if( strcasecmp(argv[nopt],"-quadrant") == 0 ){
139 smethod = "quadrant" ; nopt++ ; continue ;
140 }
141
142 if( strcasecmp(argv[nopt],"-ktaub") == 0 || strcasecmp(argv[nopt],"-taub") == 0 ){
143 smethod = "ktaub" ; nopt++ ; continue ;
144 }
145
146 if( strcasecmp(argv[nopt],"-fisher") == 0 ){ /* 12 Jan 2018 */
147 do_atanh = 1 ; nopt++ ; continue ;
148 }
149
150 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
151 prefix = argv[++nopt] ;
152 if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal value after -prefix!") ;
153 nopt++ ; continue ;
154 }
155
156 ERROR_exit("Unknown option: %s",argv[nopt]) ;
157 }
158
159 /*------------ open datasets, check for legality ------------*/
160
161 if( nopt+1 >= argc )
162 ERROR_exit("Need 2 non-option arguments on command line!?") ;
163
164 /* despite what the help says, if the 1D file is first, that's OK */
165
166 if( STRING_HAS_SUFFIX(argv[nopt],"1D") ){
167 ININFO_message("reading 1D file %s",argv[nopt]) ;
168 ysim = mri_read_1D( argv[nopt] ) ; ynam = argv[nopt] ;
169 if( ysim == NULL ) ERROR_exit("Can't read 1D file %s",argv[nopt]) ;
170 } else {
171 ININFO_message("reading dataset file %s",argv[nopt]) ;
172 xset = THD_open_dataset( argv[nopt] ) ; xnam = argv[nopt] ;
173 if( xset == NULL ) ERROR_exit("Can't open dataset %s",argv[nopt]) ;
174 }
175
176 /* read whatever type of file (3D or 1D) we don't already have */
177
178 nopt++ ;
179 if( xset != NULL ){
180 ININFO_message("reading 1D file %s",argv[nopt]) ;
181 ysim = mri_read_1D( argv[nopt] ) ; ynam = argv[nopt] ;
182 if( ysim == NULL ) ERROR_exit("Can't read 1D file %s",argv[nopt]) ;
183 } else {
184 ININFO_message("reading dataset file %s",argv[nopt]) ;
185 xset = THD_open_dataset( argv[nopt] ) ; xnam = argv[nopt] ;
186 if( xset == NULL ) ERROR_exit("Can't open dataset %s",argv[nopt]) ;
187 }
188
189 nvals = DSET_NVALS(xset) ; /* number of time points */
190 ii = (strcasecmp(smethod,"dot")==0) ? 2 : 3 ;
191 if( nvals < ii )
192 ERROR_exit("Input dataset %s length is less than ii?!",xnam,ii) ;
193
194 if( ysim->nx < nvals )
195 ERROR_exit("1D file %s has %d time points, but dataset has %d values",
196 ynam,ysim->nx,nvals) ;
197 else if( ysim->nx > nvals )
198 WARNING_message("1D file %s has %d time points, dataset has %d",
199 ynam,ysim->nx,nvals) ;
200
201 if( mri_allzero(ysim) )
202 ERROR_exit("1D file %s is all zero!",ynam) ;
203
204 if( ysim->ny > 1 )
205 INFO_message("1D file %s has %d columns: correlating with ALL of them!",
206 ynam,ysim->ny) ;
207
208 if( strcasecmp(smethod,"dot") == 0 && do_atanh ){
209 WARNING_message("'-dot' turns off '-Fisher'") ; do_atanh = 0 ;
210 }
211 if( strcasecmp(smethod,"dot") == 0 && datum == MRI_short ){
212 WARNING_message("'-dot' turns off '-short'") ; datum = MRI_float ;
213 }
214
215 cset = THD_Tcorr1D( xset, mask, nmask, ysim, smethod,
216 prefix, (datum==MRI_short) , do_atanh );
217 tross_Make_History( "3dTcorr1D" , argc,argv , cset ) ;
218
219 DSET_unload(xset) ; /* no longer needful */
220
221 /* finito */
222
223 DSET_write(cset) ;
224 INFO_message("Wrote dataset: %s\n",DSET_BRIKNAME(cset)) ;
225 exit(0) ;
226 }
227