1 #include "mrilib.h"
2
3 #undef INMASK
4 #define INMASK(i) ( mask == NULL || mask[i] != 0 )
5
6 /*-------------------------------------------------------------------*/
7
8 #undef FREEIF
9 #define FREEIF(p) do{ if((p)!=NULL){free(p);(p)=NULL;} }while(0)
10
11 static int maxlag = 0 ;
12 static float *cor = NULL ;
13 static int *ncor = NULL ;
14 static int port = 0 ;
15
16 /*....................................................................*/
17
CORMAT_init(int mlag,int pp)18 void CORMAT_init( int mlag , int pp )
19 {
20 FREEIF(cor) ; FREEIF(ncor) ; maxlag = 0 ;
21 if( mlag > 0 ){
22 cor = (float *)calloc( sizeof(float) , mlag ) ;
23 ncor = (int *) calloc( sizeof(int) , mlag ) ;
24 maxlag = mlag ;
25 port = (pp < 0 || pp > 3) ? 0 : pp ;
26 }
27 return ;
28 }
29
30 /*....................................................................*/
31
CORMAT_fetch(void)32 MRI_IMAGE * CORMAT_fetch(void)
33 {
34 MRI_IMAGE *cim ; float *car ; int ii ;
35 if( maxlag <= 0 || cor == NULL ) return NULL ;
36 cim = mri_new(maxlag,1,MRI_float) ; car = MRI_FLOAT_PTR(cim) ;
37 for( ii=0 ; ii < maxlag ; ii++ )
38 if( ncor[ii] > 0 ) car[ii] = cor[ii] / ncor[ii] ;
39 return cim ;
40 }
41
42 /*....................................................................*/
43
CORMAT_add_vector(int nv,float * vv)44 void CORMAT_add_vector( int nv , float *vv )
45 {
46 int ii,itop , ktop,kk ;
47 float ss , sq ;
48
49 if( maxlag <= 0 || nv <= 2 || vv == NULL ) return ;
50
51 switch(port){
52 default:
53 case 0: THD_const_detrend ( nv , vv , NULL ); break;
54 case 1: THD_linear_detrend ( nv , vv , NULL,NULL ); break;
55 case 2: THD_quadratic_detrend ( nv , vv , NULL,NULL,NULL ); break;
56 case 3: THD_cubic_detrend ( nv , vv ); break;
57 }
58 THD_normRMS( nv , vv ) ;
59
60 ktop = MIN(maxlag,nv-1) ;
61 for( kk=1 ; kk <= ktop ; kk++ ){
62 itop = nv-kk ;
63 for( ii=0 ; ii < itop ; ii++ ){
64 cor[kk-1] += vv[ii] * vv[ii+kk] ; ncor[kk-1]++ ;
65 }
66 }
67
68 return ;
69 }
70
71 /*-------------------------------------------------------------------*/
72
73 MRI_IMAGE * mri_cormat_vector( MRI_IMARR *imar ) ;
74 MRI_IMARR * THD_get_dset_nbhd_array( THD_3dim_dataset *dset, byte *mask,
75 int xx, int yy, int zz, MCW_cluster *nbhd ) ;
76 static void vstep_print(void) ;
77
78 float_pair estimate_arma11( int nx , float * ) ;
79
80 static int nbk=0 , *bk=NULL , ntime=0 , mlag=0 , pport=0 ; /* cheap */
81
82 /*------------------------------------------------------------------------*/
83 /* Adapted from 3dLocalSVD.c and 3dErrtsCormat.c */
84
main(int argc,char * argv[])85 int main( int argc , char *argv[] )
86 {
87 THD_3dim_dataset *inset=NULL , *outset=NULL ;
88 MCW_cluster *nbhd=NULL ;
89 byte *mask=NULL ; int mask_nx,mask_ny,mask_nz , automask=0 ;
90 char *prefix="./LocalCormat" ;
91 int iarg=1 , verb=1 , ntype=0 , kk,nx,ny,nz,nxy,nxyz,nt , xx,yy,zz, vstep ;
92 float na,nb,nc , dx,dy,dz ;
93 MRI_IMARR *imar=NULL ; MRI_IMAGE *pim=NULL ;
94 int mmlag=10 , ii,jj , do_arma=0 , nvout ;
95 MRI_IMAGE *concim=NULL ; float *concar=NULL ;
96
97 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
98 printf(
99 "Usage: 3dLocalCORMAT [options] inputdataset\n"
100 "\n"
101 "Compute the correlation matrix (in time) of the input dataset,\n"
102 "up to lag given by -maxlag. The matrix is averaged over the\n"
103 "neighborhood specified by the -nbhd option, and then the entries\n"
104 "are output at each voxel in a new dataset.\n"
105 "\n"
106 "Normally, the input to this program would be the -errts output\n"
107 "from 3dDeconvolve, or the equivalent residuals from some other\n"
108 "analysis. If you input a non-residual time series file, you at\n"
109 "least should use an appropriate -polort level for detrending!\n"
110 "\n"
111 "Options:\n"
112 " -input inputdataset\n"
113 " -prefix ppp\n"
114 " -mask mset {these 2 options are}\n"
115 " -automask {mutually exclusive.}\n"
116 " -nbhd nnn [e.g., 'SPHERE(9)' for 9 mm radius]\n"
117 " -polort ppp [default = 0, which is reasonable for -errts output]\n"
118 " -concat ccc [as in 3dDeconvolve]\n"
119 " -maxlag mmm [default = 10]\n"
120 " -ARMA [estimate ARMA(1,1) parameters into last 2 sub-bricks]\n"
121 "\n"
122 "A quick hack for my own benignant purposes -- RWCox -- June 2008\n"
123 ) ;
124 PRINT_COMPILE_DATE ; exit(0) ;
125 }
126
127 /*---- official startup ---*/
128
129 PRINT_VERSION("3dLocalCormat"); mainENTRY("3dLocalCormat main"); machdep();
130 AFNI_logger("3dLocalCormat",argc,argv); AUTHOR("Zhark the Toeplitzer");
131
132 /*---- loop over options ----*/
133
134 while( iarg < argc && argv[iarg][0] == '-' ){
135
136 #if 0
137 fprintf(stderr,"argv[%d] = %s\n",iarg,argv[iarg]) ;
138 #endif
139
140 if( strcmp(argv[iarg],"-ARMA") == 0 ){
141 do_arma = 1 ; iarg++ ; continue ;
142 }
143
144 if( strcmp(argv[iarg],"-polort") == 0 ){
145 char *cpt ;
146 if( ++iarg >= argc )
147 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
148 pport = (int)strtod(argv[iarg],&cpt) ;
149 if( *cpt != '\0' )
150 WARNING_message("Illegal non-numeric value after -polort") ;
151 if( pport > 3 ){
152 pport = 3 ; WARNING_message("-polort set to 3 == max implemented") ;
153 } else if( pport < 0 ){
154 pport = 0 ; WARNING_message("-polort set to 0 == min implemented") ;
155 }
156 iarg++ ; continue ;
157 }
158
159 if( strcmp(argv[iarg],"-input") == 0 ){
160 if( inset != NULL ) ERROR_exit("Can't have two -input options") ;
161 if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;
162 inset = THD_open_dataset( argv[iarg] ) ;
163 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
164 iarg++ ; continue ;
165 }
166
167 if( strcmp(argv[iarg],"-prefix") == 0 ){
168 if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix'") ;
169 prefix = strdup(argv[iarg]) ;
170 iarg++ ; continue ;
171 }
172
173 if( strcmp(argv[iarg],"-mask") == 0 ){
174 THD_3dim_dataset *mset ; int mmm ;
175 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
176 if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;
177 mset = THD_open_dataset( argv[iarg] ) ;
178 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
179 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
180 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
181 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
182 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;
183 mmm = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
184 INFO_message("Number of voxels in mask = %d",mmm) ;
185 if( mmm < 2 ) ERROR_exit("Mask is too small to process") ;
186 iarg++ ; continue ;
187 }
188
189 if( strcmp(argv[iarg],"-automask") == 0 ){
190 if( mask != NULL ) ERROR_exit("Can't have -automask and -mask") ;
191 automask = 1 ;
192 iarg++ ; continue ;
193 }
194
195 if( strcmp(argv[iarg],"-nbhd") == 0 ){
196 char *cpt ;
197 if( ntype > 0 ) ERROR_exit("Can't have 2 '-nbhd' options") ;
198 if( ++iarg >= argc ) ERROR_exit("Need argument after '-nbhd'") ;
199
200 cpt = argv[iarg] ;
201 if( strncasecmp(cpt,"SPHERE",6) == 0 ){
202 sscanf( cpt+7 , "%f" , &na ) ;
203 if( na == 0.0f ) ERROR_exit("Can't have a SPHERE of radius 0") ;
204 ntype = NTYPE_SPHERE ;
205 } else if( strncasecmp(cpt,"RECT",4) == 0 ){
206 sscanf( cpt+5 , "%f,%f,%f" , &na,&nb,&nc ) ;
207 if( na == 0.0f && nb == 0.0f && nc == 0.0f )
208 ERROR_exit("'RECT(0,0,0)' is not a legal neighborhood") ;
209 ntype = NTYPE_RECT ;
210 } else if( strncasecmp(cpt,"RHDD",4) == 0 ){
211 sscanf( cpt+5 , "%f" , &na ) ;
212 if( na == 0.0f ) ERROR_exit("Can't have a RHDD of radius 0") ;
213 ntype = NTYPE_RHDD ;
214 } else {
215 ERROR_exit("Unknown -nbhd shape: '%s'",cpt) ;
216 }
217 iarg++ ; continue ;
218 }
219
220 if( strcmp(argv[iarg],"-maxlag") == 0 ){
221 if( ++iarg >= argc )
222 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
223 mmlag = (int)strtod(argv[iarg],NULL) ;
224 iarg++ ; continue ;
225 }
226
227 if( strcmp(argv[iarg],"-concat") == 0 ){
228 if( concim != NULL )
229 ERROR_exit("Can't have two %s options!",argv[iarg]) ;
230 if( ++iarg >= argc )
231 ERROR_exit("Need argument after option %s",argv[iarg-1]) ;
232 concim = mri_read_1D( argv[iarg] ) ;
233 if( concim == NULL )
234 ERROR_exit("Can't read -concat file '%s'",argv[iarg]) ;
235 if( concim->nx < 2 )
236 ERROR_exit("-concat file '%s' must have at least 2 entries!",
237 argv[iarg]) ;
238 concar = MRI_FLOAT_PTR(concim) ;
239 for( ii=1 ; ii < concim->nx ; ii++ )
240 if( (int)concar[ii-1] >= (int)concar[ii] )
241 ERROR_exit("-concat file '%s' is not ordered increasingly!",
242 argv[iarg]) ;
243 iarg++ ; continue ;
244 }
245
246 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
247
248 } /*--- end of loop over options ---*/
249
250 if( do_arma && mmlag > 0 && mmlag < 5 )
251 ERROR_exit("Can't do -ARMA with -maxlag %d",mmlag) ;
252
253 /*---- deal with input dataset ----*/
254
255 if( inset == NULL ){
256 if( iarg >= argc ) ERROR_exit("No input dataset on command line?") ;
257 inset = THD_open_dataset( argv[iarg] ) ;
258 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
259 }
260 ntime = DSET_NVALS(inset) ;
261 if( ntime < 9 )
262 ERROR_exit("Must have at least 9 values per voxel") ;
263
264 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
265
266 if( mask != NULL ){
267 if( mask_nx != DSET_NX(inset) ||
268 mask_ny != DSET_NY(inset) ||
269 mask_nz != DSET_NZ(inset) )
270 ERROR_exit("-mask dataset grid doesn't match input dataset") ;
271
272 } else if( automask ){
273 int mmm ;
274 mask = THD_automask( inset ) ;
275 if( mask == NULL )
276 ERROR_message("Can't create -automask from input dataset?") ;
277 mmm = THD_countmask( DSET_NVOX(inset) , mask ) ;
278 INFO_message("Number of voxels in automask = %d",mmm) ;
279 if( mmm < 2 ) ERROR_exit("Automask is too small to process") ;
280 }
281
282 /*-- set up blocks of continuous time data --*/
283
284 if( DSET_IS_TCAT(inset) ){
285 if( concim != NULL ){
286 WARNING_message("Ignoring -concat, since dataset is auto-catenated") ;
287 mri_free(concim) ;
288 }
289 concim = mri_new(inset->tcat_num,1,MRI_float) ;
290 concar = MRI_FLOAT_PTR(concim) ;
291 concar[0] = 0.0 ;
292 for( ii=0 ; ii < inset->tcat_num-1 ; ii++ )
293 concar[ii+1] = concar[ii] + inset->tcat_len[ii] ;
294 } else if( concim == NULL ){
295 concim = mri_new(1,1,MRI_float) ;
296 concar = MRI_FLOAT_PTR(concim) ; concar[0] = 0 ;
297 }
298 nbk = concim->nx ;
299 bk = (int *)malloc(sizeof(int)*(nbk+1)) ;
300 for( ii=0 ; ii < nbk ; ii++ ) bk[ii] = (int)concar[ii] ;
301 bk[nbk] = ntime ;
302 mri_free(concim) ;
303 mlag = DSET_NVALS(inset) ;
304 for( ii=0 ; ii < nbk ; ii++ ){
305 jj = bk[ii+1]-bk[ii] ; if( jj < mlag ) mlag = jj ;
306 if( bk[ii] < 0 || jj < 9 )
307 ERROR_exit("something is rotten in the dataset run lengths") ;
308 }
309 mlag-- ;
310 if( mmlag > 0 && mlag > mmlag ) mlag = mmlag ;
311 else INFO_message("Max lag set to %d",mlag) ;
312
313 if( do_arma && mlag < 5 )
314 ERROR_exit("Can't do -ARMA with maxlag=%d",mlag) ;
315
316 /*---- create neighborhood (as a cluster) -----*/
317
318 if( ntype <= 0 ){ /* default neighborhood */
319 ntype = NTYPE_SPHERE ; na = -1.01f ;
320 INFO_message("Using default neighborhood = self + 6 neighbors") ;
321 }
322
323 switch( ntype ){
324 default:
325 ERROR_exit("WTF? ntype=%d",ntype) ;
326
327 case NTYPE_SPHERE:{
328 if( na < 0.0f ){ dx = dy = dz = 1.0f ; na = -na ; }
329 else { dx = fabsf(DSET_DX(inset)) ;
330 dy = fabsf(DSET_DY(inset)) ;
331 dz = fabsf(DSET_DZ(inset)) ; }
332 nbhd = MCW_spheremask( dx,dy,dz , na ) ;
333 }
334 break ;
335
336 case NTYPE_RECT:{
337 if( na < 0.0f ){ dx = 1.0f; na = -na; } else dx = fabsf(DSET_DX(inset));
338 if( nb < 0.0f ){ dy = 1.0f; nb = -nb; } else dy = fabsf(DSET_DY(inset));
339 if( nc < 0.0f ){ dz = 1.0f; nc = -nc; } else dz = fabsf(DSET_DZ(inset));
340 nbhd = MCW_rectmask( dx,dy,dz , na,nb,nc ) ;
341 }
342 break ;
343
344 case NTYPE_RHDD:{
345 if( na < 0.0f ){ dx = dy = dz = 1.0f ; na = -na ; }
346 else { dx = fabsf(DSET_DX(inset)) ;
347 dy = fabsf(DSET_DY(inset)) ;
348 dz = fabsf(DSET_DZ(inset)) ; }
349 nbhd = MCW_rhddmask( dx,dy,dz , na ) ;
350 }
351 break ;
352 }
353 MCW_radsort_cluster( nbhd , dx,dy,dz ) ; /* 26 Feb 2008 */
354
355 INFO_message("Neighborhood comprises %d voxels",nbhd->num_pt) ;
356
357 /** create output dataset **/
358
359 outset = EDIT_empty_copy(inset) ;
360 nvout = mlag ; if( do_arma ) nvout += 2 ;
361 EDIT_dset_items( outset,
362 ADN_prefix , prefix,
363 ADN_brick_fac, NULL ,
364 ADN_nvals , nvout ,
365 ADN_ntt , nvout ,
366 ADN_none );
367 tross_Copy_History( inset , outset ) ;
368 tross_Make_History( "3dLocalCormat" , argc,argv , outset ) ;
369 for( kk=0 ; kk < nvout ; kk++ )
370 EDIT_substitute_brick( outset , kk , MRI_float , NULL ) ;
371
372 nx = DSET_NX(outset) ;
373 ny = DSET_NY(outset) ; nxy = nx*ny ;
374 nz = DSET_NZ(outset) ; nxyz = nxy*nz ;
375 vstep = (verb && nxyz > 999) ? nxyz/50 : 0 ;
376 if( vstep ) fprintf(stderr,"++ voxel loop: ") ;
377
378 /** actually do the long long slog through all voxels **/
379
380 for( kk=0 ; kk < nxyz ; kk++ ){
381 if( vstep && kk%vstep==vstep-1 ) vstep_print() ;
382 if( !INMASK(kk) ) continue ;
383 IJK_TO_THREE( kk , xx,yy,zz , nx,nxy ) ;
384 imar = THD_get_dset_nbhd_array( inset , mask , xx,yy,zz , nbhd ) ;
385 if( imar == NULL ) continue ;
386 pim = mri_cormat_vector(imar) ; DESTROY_IMARR(imar) ;
387 if( pim == NULL ) continue ;
388 THD_insert_series( kk, outset, pim->nx, MRI_float, MRI_FLOAT_PTR(pim), 0 ) ;
389
390 if( do_arma ){ /* estimate ARMA(1,1) params and store those, too */
391 float_pair ab ;
392 float *aa=DSET_ARRAY(outset,mlag), *bb=DSET_ARRAY(outset,mlag+1) ;
393 ab = estimate_arma11( pim->nx , MRI_FLOAT_PTR(pim) ) ;
394 aa[kk] = ab.a ; bb[kk] = ab.b ;
395 }
396
397 mri_free(pim) ;
398 }
399 if( vstep ) fprintf(stderr,"\n") ;
400
401 DSET_delete(inset) ;
402 DSET_write(outset) ;
403 WROTE_DSET(outset) ;
404
405 exit(0) ;
406 }
407
408 /*------------------------------------------------------------------------*/
409
THD_get_dset_nbhd_array(THD_3dim_dataset * dset,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)410 MRI_IMARR * THD_get_dset_nbhd_array( THD_3dim_dataset *dset, byte *mask,
411 int xx, int yy, int zz, MCW_cluster *nbhd )
412 {
413 MRI_IMARR *imar ;
414 int nvox, *ivox , nx,ny,nz , nxy,nxyz , npt, aa,bb,cc,kk,ii ;
415
416 nx = DSET_NX(dset) ;
417 ny = DSET_NY(dset) ; nxy = nx*ny ;
418 nz = DSET_NZ(dset) ; nxyz = nxy*nz ; npt = nbhd->num_pt ;
419
420 kk = xx + yy*nx + zz*nxy ;
421 if( kk < 0 || kk >= nxyz || !INMASK(kk) ) return NULL ;
422
423 ivox = (int *)malloc(sizeof(int)*npt) ; nvox = 0 ;
424 for( ii=0 ; ii < npt ; ii++ ){
425 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
426 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
427 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
428 kk = aa + bb*nx + cc*nxy ;
429 if( INMASK(kk) ) ivox[nvox++] = kk ;
430 }
431
432 imar = THD_extract_many_series( nvox, ivox, dset ) ;
433 free(ivox) ; return imar ;
434 }
435
436 /*------------------------------------------------------------------------*/
437
mri_cormat_vector(MRI_IMARR * imar)438 MRI_IMAGE * mri_cormat_vector( MRI_IMARR *imar )
439 {
440 int nx , nvec , ii,jj ;
441 double *amat , *umat , *vmat , *sval ;
442 float *far ; MRI_IMAGE *tim ;
443
444 if( imar == NULL ) return NULL ;
445 nvec = IMARR_COUNT(imar) ; if( nvec < 1 ) return NULL ;
446 nx = IMARR_SUBIM(imar,0)->nx ; if( nx < 9 ) return NULL ;
447
448 CORMAT_init( mlag , pport ) ;
449
450 for( jj=0 ; jj < nvec ; jj++ ){
451 tim = IMARR_SUBIM(imar,jj) ;
452 far = MRI_FLOAT_PTR(tim) ;
453 for( ii=0 ; ii < nbk ; ii++ )
454 CORMAT_add_vector( bk[ii+1]-bk[ii] , far + bk[ii] ) ;
455 }
456
457 tim = CORMAT_fetch() ; return tim ;
458 }
459
460 /*------------------------------------------------------------------------*/
461
vstep_print(void)462 static void vstep_print(void)
463 {
464 static char xx[10] = "0123456789" ; static int vn=0 ;
465 fprintf(stderr , "%c" , xx[vn%10] ) ;
466 if( vn%10 == 9) fprintf(stderr,".") ;
467 vn++ ;
468 }
469
470 /*------------------------------------------------------------------------*/
471 /* Estimation of ARMA(1,1) coefficients a and b from correlations. */
472
473 static int nc ; /* number of correlations */
474 static float *rc ; /* correlation vector */
475
arma11func(int np,double * par)476 double arma11func( int np , double *par ) /* cost = weighted least squares */
477 {
478 register double a,b,sum,rr ; register int kk ;
479
480 a = par[0] ; b = par[1] ;
481 rr = (b+a)*(1.0+a*b)/(1.0+b*(b+2.0*a)) ; /* correlation at lag=1 */
482 sum = (rr-rc[0])*(rr-rc[0]) * (0.2+rc[0]*rc[0]) ;
483 for( kk=1 ; kk < nc ; kk++ ){
484 rr *= a ; /* correlation at lag=kk+1 */
485 sum += (rr-rc[kk])*(rr-rc[kk]) * (0.2+rc[kk]*rc[kk]) ;
486 }
487 return sum ;
488 }
489
estimate_arma11(int nx,float * cc)490 float_pair estimate_arma11( int nx , float *cc )
491 {
492 double x[2] , xbot[2],xtop[2] ;
493 float_pair ab ;
494
495 nc = nx ; rc = cc ; /* set static variables for arma11func() */
496 x[0] = x[1] = 0.1 ; /* initial guesses for a and b */
497 xbot[0] = 0.0 ; xtop[0] = 0.8 ; /* range of allowed values for a and b */
498 xbot[1] = 0.0 ; xtop[1] = 0.8 ;
499
500 (void)powell_newuoa_constrained( 2 , x , NULL , xbot , xtop ,
501 13 , 5 , 2 ,
502 0.10 , 0.005 , 222 , arma11func ) ;
503 ab.a = (float)x[0] ;
504 ab.b = (float)x[1] ; return ab ;
505 }
506