1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 #define SPEARMAN 1
8 #define QUADRANT 2
9 #define PEARSON 3
10 #define KTAUB 4
11 #define DOT 5
12
13 #undef MYatanh
14 #define MYatanh(x) ( ((x)<-0.999329f) ? -4.0f \
15 :((x)>+0.999329f) ? +4.0f : atanhf(x) )
16
17 /*-------------------------------------------------------------*/
THD_dotprod(int nx,float * xx,float * yy)18 float THD_dotprod( int nx , float *xx , float *yy ){
19 int jj ; float sum=0.0f ;
20 for( jj=0 ; jj < nx ; jj++ ) sum += xx[jj]*yy[jj] ;
21 return sum ;
22 }
23 /*-------------------------------------------------------------*/
24
THD_Tcorr1D(THD_3dim_dataset * xset,byte * mask,int nmask,MRI_IMAGE * ysim,char * smethod,char * prefix,int do_short,int do_atanh)25 THD_3dim_dataset *THD_Tcorr1D(THD_3dim_dataset *xset, byte *mask, int nmask,
26 MRI_IMAGE *ysim,
27 char *smethod, char *prefix, int do_short , int do_atanh )
28 {
29 THD_3dim_dataset *cset = NULL;
30 int method=PEARSON ;
31 int ny, kk, datum=MRI_float ; char fmt[32] ; float cfac=0.0f,sfac=0.0f ;
32 float (*corfun)(int,float *,float *) = NULL ; /* ptr to corr function */
33 int nvox , nvals , ii;
34 int nconst=0 ;
35 char blab[128], alab[192] ; /* 31 Aug 2021 */
36
37 ENTRY("THD_Tcorr1D");
38
39 if( do_short ) datum = MRI_short ; /* 30 Jan 2017 */
40
41 if (!smethod || smethod[0] == '\0') {
42 method = PEARSON;
43 } else if (!strcasecmp(smethod,"pearson")) {
44 method = PEARSON;
45 } else if (!strcasecmp(smethod,"spearman")) {
46 method = SPEARMAN;
47 } else if (!strcasecmp(smethod,"quadrant")) {
48 method = QUADRANT;
49 } else if (!strcasecmp(smethod,"ktaub")) {
50 method = KTAUB;
51 } else if (!strcasecmp(smethod,"dot")) {
52 method = DOT;
53 if( do_atanh ){
54 WARNING_message("Method %s disables the Fishter transformation",smethod) ;
55 do_atanh = 0 ;
56 }
57 } else {
58 WARNING_message("Bad value %s for correlation method ==> PEARSON", smethod);
59 method = PEARSON;
60 }
61
62 if (!prefix) prefix = "Tcorr1D";
63
64 nvals = DSET_NVALS(xset) ; /* number of time points */
65
66 ii = (method==DOT) ? 2 : 3 ;
67 if( nvals < ii )
68 ERROR_exit("Input dataset length (%d) is less than %d?!", nvals,ii) ;
69
70 if( ysim->nx < nvals )
71 ERROR_exit("ysim has %d time points, but dataset has %d values",
72 ysim->nx,nvals) ;
73 else if( ysim->nx > nvals )
74 WARNING_message("ysim has %d time points, dataset has %d",
75 ysim->nx,nvals) ;
76
77 if( mri_allzero(ysim) )
78 ERROR_exit("ysim is all zero!") ;
79
80 ny = ysim->ny ;
81 if( ny > 1 )
82 INFO_message("ysim has %d columns: correlating with ALL of them!",
83 ny) ;
84
85 ININFO_message("loading dataset %s into memory",DSET_BRIKNAME(xset)) ;
86 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
87
88 nvox = DSET_NVOX(xset) ;
89 if( mask == NULL ) nmask = nvox ;
90
91 /*-- create output dataset --*/
92
93 cset = EDIT_empty_copy( xset ) ;
94 EDIT_dset_items( cset ,
95 ADN_prefix , prefix ,
96 ADN_nvals , ny ,
97 ADN_ntt , 0 , /* no time axis */
98 ADN_type , HEAD_FUNC_TYPE ,
99 ADN_func_type , FUNC_BUCK_TYPE ,
100 ADN_none ) ;
101
102 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(cset)) )
103 ERROR_exit("Output dataset %s already exists!",DSET_HEADNAME(cset)) ;
104
105 if( ny <= 10 ) kk = 1 ; /* number of digits for */
106 else if( ny <= 100 ) kk = 2 ; /* brick label string */
107 else if( ny <= 1000 ) kk = 3 ;
108 else if( ny <= 10000 ) kk = 4 ;
109 else kk = 5 ; /* unlikely! */
110 switch( method ){ /* make brick label string format */
111 default:
112 case PEARSON: sprintf(fmt,"PearCorr#%%0%dd",kk) ; break ;
113 case SPEARMAN: sprintf(fmt,"SpmnCorr#%%0%dd",kk) ; break ;
114 case QUADRANT: sprintf(fmt,"QuadCorr#%%0%dd",kk) ; break ;
115 case KTAUB: sprintf(fmt,"TaubCorr#%%0%dd",kk) ; break ;
116 case DOT: sprintf(fmt, "DotProd#%%0%dd",kk) ; break ;
117 }
118 if( datum == MRI_short ){
119 cfac = (do_atanh) ? 0.000125f : 0.0001f ; /* scale factor for -short */
120 sfac = 1.0f/cfac + 0.111f ;
121 }
122
123 /* for each sub-brick in output file:
124 create empty volume
125 mark volume with statistic code
126 set volume scale factor
127 set volume label */
128
129 for( kk=0 ; kk < ny ; kk++ ){
130 EDIT_substitute_brick(cset,kk,datum,NULL) ; /* make brick */
131 if( do_atanh ){
132 EDIT_BRICK_TO_FIZT(cset,kk) ;
133 } else if( method != DOT ){
134 EDIT_BRICK_TO_FICO(cset,kk,nvals,1,1) ; /* stat params */
135 } else {
136 EDIT_BRICK_TO_NOSTAT(cset,kk) ;
137 }
138 EDIT_BRICK_FACTOR(cset,kk,cfac) ; /* set brick factor */
139 sprintf(blab,fmt,kk) ; /* manufacture label */
140 if( do_atanh ){
141 sprintf(alab,"atanh_%s",blab) ;
142 } else {
143 strcpy(alab,blab) ;
144 }
145 EDIT_BRICK_LABEL(cset,kk,alab) ; /* labelize brick */
146 }
147
148 switch( method ){ /* set correlation function */
149 default:
150 case PEARSON: corfun = THD_pearson_corr ; break ;
151 case SPEARMAN: corfun = THD_spearman_corr ; break ;
152 case QUADRANT: corfun = THD_quadrant_corr ; break ;
153 case KTAUB: corfun = THD_ktaub_corr ; break ;
154 case DOT: corfun = THD_dotprod ; break ;
155 }
156
157 /* 27 Jun 2010: OpenMP-ize over columns in ysim */
158
159 AFNI_OMP_START ;
160 #pragma omp parallel if( ny > 1 )
161 { float *ysar, *xsar, *fcar = NULL, *ydar, val ;
162 int ii, kk, jj ; short *scar = NULL;
163
164 #ifdef USE_OMP
165 if( omp_get_thread_num() == 0 )
166 INFO_message("Start correlations: %d voxels X %d time series(%d); %d threads",
167 nmask , ny , nvals , omp_get_num_threads() ) ;
168 #else
169 INFO_message("Start correlations: %d voxels X %d time series(%d)",nmask,ny,nvals) ;
170 #endif
171
172 ydar = (float *)malloc(sizeof(float)*nvals) ; /* 1D data duplicate */
173 xsar = (float *)malloc(sizeof(float)*nvals) ; /* 3D data duplicate */
174
175 /* 27 Jun 2010: loop over columns in ysim */
176
177 #pragma omp for
178 for( kk=0 ; kk < ny ; kk++ ){ /* loop over ysim columns */
179 if( datum == MRI_short ) scar = DSET_ARRAY(cset,kk) ; /* output array */
180 else fcar = DSET_ARRAY(cset,kk) ;
181 ysar = MRI_FLOAT_PTR(ysim) + (kk * ysim->nx) ; /* 1D data pointer */
182
183 /* loop over voxels, correlate */
184
185 for( ii=0 ; ii < nvox ; ii++ ){
186
187 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip this'n */
188
189 /* get time series to correlate */
190
191 (void)THD_extract_array(ii,xset,0,xsar) ; /* 3D data */
192 for( jj=1 ; jj < nvals && xsar[jj]==xsar[0] ; jj++ ) ; /* nada */
193 if( jj == nvals ){ /* data was constant */
194 #pragma omp atomic
195 nconst++ ;
196 continue ;
197 }
198 for( jj=0 ; jj < nvals ; jj++ ) ydar[jj] = ysar[jj] ; /* 1D data */
199
200 val = corfun( nvals , xsar , ydar ) ; /* !! correlate !! */
201 if( do_atanh ) val = MYatanh(val) ;
202
203 if( datum == MRI_short ) scar[ii] = (short)(sfac*val) ;
204 else fcar[ii] = val ;
205
206 } /* end of loop over voxels */
207
208 #pragma omp critical
209 { if( ny > 1 ) fprintf(stderr,"[%d]",kk) ; }
210 } /* end of loop over ysim columns */
211
212 free(ydar) ; free(xsar) ;
213 } /* end OpenMP */
214 AFNI_OMP_END ;
215
216 if( ny > 1 ){ fprintf(stderr,"\n") ; nconst /= ny ; }
217 if( nconst > 0 )
218 WARNING_message("THD_Tcorr1D: %d voxel%s skipped because were constant in time",
219 nconst , (nconst==1) ? "\0" : "s" ) ;
220
221 RETURN(cset);
222 }
223