1 #include "mrilib.h"
2
3 /*-------------------------------------------------------------------*/
4
5 #undef FREEIF
6 #define FREEIF(p) do{ if((p)!=NULL){free(p);(p)=NULL;} }while(0)
7
8 static int maxlag = 0 ;
9 static float *cor = NULL ;
10 static int *ncor = NULL ;
11 static int port = 0 ;
12
13 /*....................................................................*/
14
CORMAT_init(int mlag,int pp)15 void CORMAT_init( int mlag , int pp )
16 {
17 FREEIF(cor) ; FREEIF(ncor) ; maxlag = 0 ;
18 if( mlag > 0 ){
19 cor = (float *)calloc( sizeof(float) , mlag ) ;
20 ncor = (int *) calloc( sizeof(int) , mlag ) ;
21 maxlag = mlag ;
22 port = (pp < 0 || pp > 3) ? 0 : pp ;
23 }
24 return ;
25 }
26
27 /*....................................................................*/
28
CORMAT_fetch(void)29 MRI_IMAGE * CORMAT_fetch(void)
30 {
31 MRI_IMAGE *cim ; float *car ; int ii ;
32 if( maxlag <= 0 || cor == NULL ) return NULL ;
33 cim = mri_new(maxlag,1,MRI_float) ; car = MRI_FLOAT_PTR(cim) ;
34 for( ii=0 ; ii < maxlag ; ii++ )
35 if( ncor[ii] > 0 ) car[ii] = cor[ii] / ncor[ii] ;
36 return cim ;
37 }
38
39 /*....................................................................*/
40
CORMAT_add_vector(int nv,float * vv)41 void CORMAT_add_vector( int nv , float *vv )
42 {
43 int ii,itop , ktop,kk ;
44 float ss , sq ;
45
46 if( maxlag <= 0 || nv <= 2 || vv == NULL ) return ;
47
48 switch(port){
49 default:
50 case 0: THD_const_detrend ( nv , vv , NULL ); break;
51 case 1: THD_linear_detrend ( nv , vv , NULL,NULL ); break;
52 case 2: THD_quadratic_detrend ( nv , vv , NULL,NULL,NULL ); break;
53 case 3: THD_cubic_detrend ( nv , vv ); break;
54 }
55 THD_normRMS( nv , vv ) ;
56
57 ktop = MIN(maxlag,nv-1) ;
58 for( kk=1 ; kk <= ktop ; kk++ ){
59 itop = nv-kk ;
60 for( ii=0 ; ii < itop ; ii++ ){
61 cor[kk-1] += vv[ii] * vv[ii+kk] ; ncor[kk-1]++ ;
62 }
63 }
64
65 return ;
66 }
67
68 /*-------------------------------------------------------------------*/
69
Syntax(void)70 void Syntax(void)
71 {
72 printf(
73 "Usage: 3dErrtsCormat [options] dset\n"
74 "\n"
75 "Computes the correlation (not covariance) matrix corresponding\n"
76 "to the residual (or error) time series in 'dset', which will\n"
77 "usually be the '-errts' output from 3dDeconvolve. The output\n"
78 "is a 1D file of the Toeplitz entries (to stdout).\n"
79 "\n"
80 "Options:\n"
81 " -concat rname = as in 3dDeconvolve\n"
82 " -input dset = alternate way of telling what dataset to read\n"
83 " -mask mset = mask dataset\n"
84 " -maxlag mm = set maximum lag\n"
85 " -polort pp = set polort level (default=0)\n"
86 "\n"
87 "-- RWCox -- June 2008 -- for my own pleasant purposes\n"
88 "-- Also see program 3dLocalCormat to do this on each voxel,\n"
89 " and optionally estimate the ARMA(1,1) model parameters.\n"
90 ) ;
91 PRINT_COMPILE_DATE ; exit(0) ;
92 }
93
94 /*-------------------------------------------------------------------*/
95
main(int argc,char * argv[])96 int main( int argc , char *argv[] )
97 {
98 MRI_IMAGE *concim=NULL, *corim ;
99 float *concar=NULL , *corar ;
100 THD_3dim_dataset *inset=NULL , *mset=NULL ;
101 int iarg , ii,jj,nvox,ntime , nbk,*bk,mlag , mmlag=0,pport=0 ;
102 byte *mmm ; float *tar ;
103
104 /*-- basic startup stuff --*/
105
106 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Syntax() ;
107
108 mainENTRY("3dErrtsCormat") ; machdep() ;
109 PRINT_VERSION("3dErrtsCormat") ; AUTHOR("Zhark the Correlator") ;
110 AFNI_logger("3dErrtsCormat",argc,argv) ;
111
112 /*-- process command line options --*/
113
114 iarg = 1 ;
115 while( iarg < argc && argv[iarg][0] == '-' ){
116
117 if( strcmp(argv[iarg],"-polort") == 0 ){
118 char *qpt ;
119 if( ++iarg >= argc )
120 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
121 pport = (int)strtod(argv[iarg],&qpt) ;
122 if( *qpt != '\0' ) WARNING_message("Illegal non-numeric value after -polort") ;
123 iarg++ ; continue ;
124 }
125
126 if( strcmp(argv[iarg],"-maxlag") == 0 ){
127 if( ++iarg >= argc )
128 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
129 mmlag = (int)strtod(argv[iarg],NULL) ;
130 iarg++ ; continue ;
131 }
132
133 if( strcmp(argv[iarg],"-concat") == 0 ){
134 if( concim != NULL )
135 ERROR_exit("Can't have two %s options!",argv[iarg]) ;
136 if( ++iarg >= argc )
137 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
138 concim = mri_read_1D( argv[iarg] ) ;
139 if( concim == NULL )
140 ERROR_exit("Can't read -concat file '%s'",argv[iarg]) ;
141 if( concim->nx < 2 )
142 ERROR_exit("-concat file '%s' must have at least 2 entries!",
143 argv[iarg]) ;
144 concar = MRI_FLOAT_PTR(concim) ;
145 for( ii=1 ; ii < concim->nx ; ii++ )
146 if( (int)concar[ii-1] >= (int)concar[ii] )
147 ERROR_exit("-concat file '%s' is not ordered increasingly!",
148 argv[iarg]) ;
149 iarg++ ; continue ;
150 }
151
152 if( strcmp(argv[iarg],"-input") == 0 ){
153 if( inset != NULL )
154 ERROR_exit("Can't have two %s options!",argv[iarg]) ;
155 if( ++iarg >= argc )
156 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
157 inset = THD_open_dataset( argv[iarg] ) ;
158 if( inset == NULL )
159 ERROR_exit("Can't open -input dataset '%s'",argv[iarg]) ;
160 if( DSET_NVALS(inset) < 9 )
161 ERROR_exit("Dataset '%s' has only %d time points",
162 argv[iarg],DSET_NVALS(inset) ) ;
163 iarg++ ; continue ;
164 }
165
166 if( strcmp(argv[iarg],"-mask") == 0 ){
167 if( mset != NULL )
168 ERROR_exit("Can't have two %s options!",argv[iarg]) ;
169 if( ++iarg >= argc )
170 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
171 mset = THD_open_dataset( argv[iarg] ) ;
172 if( mset == NULL )
173 ERROR_exit("Can't open -input dataset '%s'",argv[iarg]) ;
174 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
175 iarg++ ; continue ;
176 }
177
178 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
179 }
180
181 /*-- get input dataset --*/
182
183 if( iarg == argc && inset == NULL )
184 ERROR_exit("No input dataset on command line?") ;
185
186 if( inset == NULL ){
187 inset = THD_open_dataset( argv[iarg] ) ;
188 if( inset == NULL )
189 ERROR_exit("Can't open -input dataset '%s'",argv[iarg]) ;
190 if( DSET_NVALS(inset) < 9 )
191 ERROR_exit("Dataset '%s' has only %d time points",
192 argv[iarg],DSET_NVALS(inset) ) ;
193 }
194 DSET_load(inset); CHECK_LOAD_ERROR(inset);
195
196 /*-- create byte mask array --*/
197
198 nvox = DSET_NVOX(inset) ; ntime = DSET_NVALS(inset) ;
199
200 if( mset != NULL ){
201 if( !EQUIV_GRIDS(inset,mset) )
202 ERROR_exit("Input and mask datasets don't have the same 3D grid!") ;
203 mmm = THD_makemask( mset , 0 , 1.0f,-1.0f ) ;
204 INFO_message("%d voxels in mask (out of %d total)",
205 THD_countmask(nvox,mmm),nvox) ;
206 DSET_unload( mset ) ;
207 } else {
208 mmm = (byte *)malloc(sizeof(byte)*nvox) ;
209 memset( mmm , 1 , sizeof(byte)*nvox ) ;
210 INFO_message("Using all %d voxels in the dataset",nvox) ;
211 }
212
213 /*-- set up blocks of continuous time data --*/
214
215 if( DSET_IS_TCAT(inset) ){
216 if( concim != NULL ){
217 WARNING_message("Ignoring -concat, since dataset is auto-catenated") ;
218 mri_free(concim) ;
219 }
220 concim = mri_new(inset->tcat_num,1,MRI_float) ;
221 concar = MRI_FLOAT_PTR(concim) ;
222 concar[0] = 0.0 ;
223 for( ii=0 ; ii < inset->tcat_num-1 ; ii++ )
224 concar[ii+1] = concar[ii] + inset->tcat_len[ii] ;
225 } else if( concim == NULL ){
226 concim = mri_new(1,1,MRI_float) ;
227 concar = MRI_FLOAT_PTR(concim) ; concar[0] = 0 ;
228 }
229 nbk = concim->nx ;
230 bk = (int *)malloc(sizeof(int)*(nbk+1)) ;
231 for( ii=0 ; ii < nbk ; ii++ ) bk[ii] = (int)concar[ii] ;
232 bk[nbk] = ntime ;
233 mri_free(concim) ;
234 mlag = DSET_NVALS(inset) ;
235 for( ii=0 ; ii < nbk ; ii++ ){
236 jj = bk[ii+1]-bk[ii] ; if( jj < mlag ) mlag = jj ;
237 if( bk[ii] < 0 || jj < 9 )
238 ERROR_exit("something is rotten in the dataset run lengths") ;
239 }
240 mlag-- ;
241 if( mmlag > 0 && mlag > mmlag ) mlag = mmlag ;
242 else INFO_message("Max lag set to %d",mlag) ;
243
244 /*-- do some work --*/
245
246 CORMAT_init(mlag,pport) ;
247 tar = (float *)malloc(sizeof(float)*ntime) ;
248
249 for( jj=0 ; jj < nvox ; jj++ ){
250 if( !mmm[jj] ) continue ;
251 ii = THD_extract_array( jj , inset , 0 , tar ) ;
252 if( ii < 0 ){ ERROR_message("bad data at voxel index #%d?",jj); continue; }
253 for( ii=0 ; ii < nbk ; ii++ )
254 CORMAT_add_vector( bk[ii+1]-bk[ii] , tar + bk[ii] ) ;
255 }
256
257 /*-- finish up --*/
258
259 DSET_unload(inset) ;
260 corim = CORMAT_fetch() ;
261 if( corim == NULL ) ERROR_exit("CORMAT_fetch failed!?") ;
262 corar = MRI_FLOAT_PTR(corim) ;
263 printf(" 1.0\n") ;
264 for( ii=0 ; ii < corim->nx ; ii++ ) printf(" %.5f\n",corar[ii]) ;
265
266 exit(0) ;
267 }
268