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