1 #include "mrilib.h"
2 
3 #undef  EM
4 #define EM(s) ERROR_message("InstaCorr setup bad: %s",(s))
5 
6 /*---------------------------------------------------------------------------*/
7 /* Return 1 if methods uses timeseries normalization, 0 otherwise */
8 
THD_instacorr_cmeth_needs_norm(int cmeth)9 int THD_instacorr_cmeth_needs_norm(int cmeth) {
10    switch (cmeth) {
11       case NBISTAT_EUCLIDIAN_DIST:
12       case NBISTAT_CITYBLOCK_DIST:
13          return(0);
14       default:
15          return(1);
16    }
17    return(1);
18 }
19 
20 /*---------------------------------------------------------------------------*/
21 /* Read and pre-process time series for InstaCorr [moved here 10 Nov 2011].  */
22 
THD_instacorr_tsprep(ICOR_setup * iset,THD_3dim_dataset * dset)23 static MRI_vectim * THD_instacorr_tsprep( ICOR_setup *iset , THD_3dim_dataset *dset )
24 {
25    MRI_vectim *mv ;
26    int iv , nmmm , ntime ;
27    float **dvec , **gvec=NULL ; int ngvec=0 ;
28    int ngext=0 ; float *uvec=NULL ;
29 
30 ENTRY("THD_instacorr_tsprep") ;
31 
32    if( iset == NULL || dset == NULL ) RETURN(NULL) ;
33 
34    /*--- Extract time series for analysis ---*/
35 
36    ININFO_message("Extracting dataset time series: %s",DSET_BRIKNAME(dset)) ;
37 
38    mv = THD_dset_to_vectim_stend( dset , iset->mmm , iset->start , iset->end ) ;
39    if( mv == NULL ){
40      EM("Can't extract dataset time series!") ; RETURN(NULL) ;
41    }
42    nmmm  = mv->nvec ;
43    ntime = mv->nvals ;  /* #dataset time points to process */
44 
45    if( iset->despike ){
46      int_pair ip ;
47      ININFO_message("- Testing time series for spikes") ;
48      ip = THD_vectim_despike9( mv ) ;
49      if( ip.i > 0 )
50        ININFO_message(" -- Removed %d spikes from %d time series",ip.j,ip.i) ;
51      else
52        ININFO_message(" -- No spikes were found") ;
53    }
54 
55    /*--- Filter time series ---*/
56 
57    ININFO_message("- Filtering %d dataset time series",nmmm) ;
58 
59    if( iset->fbot <  0.0f       ) iset->fbot = 0.0f ;
60    if( iset->ftop <= iset->fbot ) iset->ftop = 999999.9f ;  /* infinity */
61 
62    dvec = (float **)malloc(sizeof(float *)*nmmm) ;
63    for( iv=0 ; iv < nmmm ; iv++ ) dvec[iv] = VECTIM_PTR(mv,iv) ;
64 
65    if( iset->gortim != NULL ){         /** set up the extra orts **/
66      if( iset->gortim->nx < ntime ){
67 
68        /* too short to be useful */
69 
70        ERROR_message("Global Ort time series length=%d is too short!",iset->gortim->nx);
71        ERROR_message(" ==> ignoring Global Ort file") ;
72 
73      } else if( iset->gortim->nx >= ntime && iset->gortim->nx < ntime+iset->start ){
74 
75        /* too short to ignore initial points */
76 
77        if( iset->start > 0 )
78          ININFO_message("Global Ort time series length=%d: not ignoring initial points",
79                         iset->gortim->nx ) ;
80        ngvec = iset->gortim->ny ;
81        gvec  = (float **)malloc(sizeof(float *)*ngvec) ;
82        for( iv=0 ; iv < ngvec ; iv++ )
83          gvec[iv] = MRI_FLOAT_PTR(iset->gortim) + iv*iset->gortim->nx ;
84 
85      } else {
86 
87        /* long enough to ignore initial points */
88 
89        if( iset->start > 0 )
90          ININFO_message("Global Ort time series length=%d: ignoring first %d points",
91                         iset->gortim->nx , iset->start ) ;
92        ngvec = iset->gortim->ny ;
93        gvec  = (float **)malloc(sizeof(float *)*ngvec) ;
94        for( iv=0 ; iv < ngvec ; iv++ )
95          gvec[iv] = MRI_FLOAT_PTR(iset->gortim) + iv*iset->gortim->nx + iset->start ;
96 
97      }
98    }
99 
100    ngext = iset->gortnpc ;
101    if( ngext > 0 ){
102      MRI_IMAGE *qim ; MRI_IMARR *qimar ; float *qar , *svec ;
103      INIT_IMARR(qimar) ;
104      for( iv=0 ; iv < nmmm ; iv++ ){
105        qim = mri_new(ntime,1,MRI_float) ;
106        qar = MRI_FLOAT_PTR(qim) ;
107        memcpy(qar,dvec[iv],sizeof(float)*ntime) ;
108        if( iset->polort >= 0 ) THD_cubic_detrend(ntime,qar) ;
109        ADDTO_IMARR(qimar,qim) ;
110      }
111      uvec = (float *)malloc(sizeof(float)*ntime*ngext) ;
112      svec = (float *)malloc(sizeof(float)*ngext) ;
113      iv = mri_principal_vectors( qimar , ngext , svec , uvec ) ;
114      DESTROY_IMARR(qimar) ;
115      if( iv < ngext ){
116        WARNING_message("Could not compute PC vectors -- skipping :-(((") ;
117        free(uvec) ;
118      } else {
119        gvec = (float **)realloc(gvec,sizeof(float *)*(ngvec+ngext)) ;
120        for( iv=0 ; iv < ngext ; iv++ )
121          gvec[iv+ngvec] = uvec + iv*ntime ;
122        fprintf(stderr," + singular values:") ;
123        for( iv=0 ; iv < ngext ; iv++ ) fprintf(stderr," %.6g",svec[iv]) ;
124        fprintf(stderr,"\n") ;
125        free(svec) ;
126      }
127    }
128 
129    /* despite the 'bandpass' name, this is the filtering step */
130 
131    (void)THD_bandpass_OK( ntime , mv->dt , iset->fbot,iset->ftop , 1 ) ;
132 
133    iset->ndet = THD_bandpass_vectors( ntime, nmmm, dvec, mv->dt,
134                                       iset->fbot, iset->ftop, iset->polort,
135                                       ngvec+ngext, gvec ) ;
136 
137 /** ININFO_message("Filtering removed %d DOF",iset->ndet) ; **/
138 
139    free(dvec) ;
140    if( gvec != NULL ) free(gvec) ;
141    if( uvec != NULL ) free(uvec) ;
142 
143    /*--- Blur time series ---*/
144 
145    if( iset->blur > 0.0f ){
146      int nrep ; float fx,fy,fz ;
147      ININFO_message("- Starting %.1f mm blur %s" ,
148                     iset->blur , (iset->mmm==NULL) ? "\0" : "(in mask)" ) ;
149      mri_blur3D_getfac( iset->blur ,
150                         mv->dx , mv->dy , mv->dz ,
151                         &nrep , &fx , &fy , &fz ) ;
152      ININFO_message("- Spatially blurring %d dataset volumes: %d iterations",mv->nvals,nrep) ;
153      mri_blur3D_vectim( mv , iset->blur ) ;
154    }
155 
156    /*-- normalize --*/
157 
158    if( THD_instacorr_cmeth_needs_norm(iset->cmeth) ){
159      ININFO_message("- Normalizing dataset time series") ;
160      THD_vectim_normalize( mv ) ;
161    }
162 
163    RETURN(mv) ;
164 }
165 
166 /*---------------------------------------------------------------------------*/
167 /*! Return value is 0 if an error occurs, or number of voxels prepared. */
168 
THD_instacorr_prepare(ICOR_setup * iset)169 int THD_instacorr_prepare( ICOR_setup *iset )
170 {
171    int nmmm=0 ; byte *mmm=NULL ;
172 
173 ENTRY("THD_instacorr_prepare") ;
174 
175    /*-- check inputs for a reasonable amount of usefulness --*/
176 
177    if( iset == NULL || !ISVALID_DSET(iset->dset) ){ EM("no dataset"); RETURN(0); }
178    if( DSET_NVALS(iset->dset) < 4 )               { EM("nvals < 4") ; RETURN(0); }
179    if( ISVALID_DSET(iset->eset) &&
180        ( DSET_NVOX (iset->eset) != DSET_NVOX (iset->dset) ||
181          DSET_NVALS(iset->eset) != DSET_NVALS(iset->dset)   ) ){ EM("eset <> dset"); RETURN(0); }
182 
183    /*-- create mask --*/
184 
185    if( iset->mmm != NULL ){ free(iset->mmm) ; iset->mmm = NULL ; }
186 
187    INFO_message("InstaCorr preparations:") ;
188 
189    /*-- automask? --*/
190 
191    if( iset->automask ){
192      ININFO_message("Computing automask") ;
193      mmm = THD_automask(iset->dset) ;
194      if( mmm == NULL ){
195        ERROR_message("Can't create automask from '%s'?!",DSET_BRIKNAME(iset->dset)) ;
196        RETURN(0) ;
197      }
198      nmmm = THD_countmask( DSET_NVOX(iset->dset) , mmm ) ;
199      if( nmmm < 9 ){
200        ERROR_message("Automask from '%s' has %d voxels!" ,
201                      DSET_BRIKNAME(iset->dset) , nmmm ) ;
202        free(mmm) ; RETURN(0) ;
203      }
204      ININFO_message("Automask from '%s' has %d voxels",DSET_BRIKNAME(iset->dset),nmmm) ;
205      iset->mmm = mmm ;
206    }
207 
208    /*-- or read a mask dataset? ---*/
209 
210    if( iset->mmm == NULL && iset->mset != NULL ){
211      if( DSET_NVOX(iset->mset) != DSET_NVOX(iset->dset) ){
212        ERROR_message("Mask dataset '%s' doesn't match input dataset '%s'",
213                      DSET_BRIKNAME(iset->mset) , DSET_BRIKNAME(iset->dset) ) ;
214        RETURN(0) ;
215      }
216      ININFO_message("Extracting mask from mask dataset sub-brick #%d",iset->mindex) ;
217      mmm = THD_makemask( iset->mset , iset->mindex , 1.0f,-1.0f ) ;
218      if( mmm == NULL ){
219        ERROR_message("Can't make mask from '%s[%d]'" ,
220                      DSET_BRIKNAME(iset->mset) , iset->mindex ) ;
221        RETURN(0) ;
222      }
223      nmmm = THD_countmask( DSET_NVOX(iset->mset) , mmm ) ;
224      if( nmmm < 9 ){
225        ERROR_message("Mask from '%s[%d]' has %d voxels!" ,
226                      DSET_BRIKNAME(iset->mset) , iset->mindex , nmmm ) ;
227        free(mmm) ; RETURN(0) ;
228      }
229      ININFO_message("Mask from '%s[%d]' has %d voxels!" ,
230                    DSET_BRIKNAME(iset->mset) , iset->mindex , nmmm ) ;
231      iset->mmm = mmm ;
232    }
233 
234    if( iset->mmm == NULL ){
235      ININFO_message("No mask for InstaCorr") ;
236      nmmm = DSET_NVOX(iset->dset) ;
237    }
238 
239    /*--- time series preparation ---*/
240 
241    iset->mv = THD_instacorr_tsprep( iset , iset->dset ) ;
242    if( iset->mv == NULL ) RETURN(0) ;
243    iset->ev = THD_instacorr_tsprep( iset , iset->eset ) ;
244 
245    RETURN(nmmm) ;
246 }
247 
248 /*---------------------------------------------------------------------------*/
249 /* Get the seed timeseries [moved to a function 07 Oct 2014] */
250 
THD_instacorr_getseed(ICOR_setup * iset,int ijk)251 float * THD_instacorr_getseed( ICOR_setup *iset , int ijk )
252 {
253    float *tsar ; int kk ; float sblur ;
254 
255 ENTRY("THD_instacorr_getseed") ;
256 
257    if( iset == NULL || iset->mv == NULL || ijk < 0 ) RETURN(NULL) ;
258 
259    /** extract reference time series **/
260 
261    tsar = (float *)malloc(sizeof(float)*(iset->mv->nvals+iset->start)) ;
262    kk   = THD_vectim_ifind( ijk , iset->mv ) ;
263 
264    if( kk >= 0 ){ /* direct from vectim, if available */
265 
266      memcpy( tsar , VECTIM_PTR(iset->mv,kk) , sizeof(float)*iset->mv->nvals ) ;
267 
268    } else {       /* non-vectim voxel in dataset ==> must be processed */
269 
270      free(tsar) ; RETURN(NULL) ;  /* don't allow this */
271 
272    }
273 
274    /** blur the ref time series, if ordered [15 May 2009] **/
275 
276    sblur = iset->sblur ;
277 
278    if( sblur != 0.0f ){
279      int gblur = AFNI_yesenv("AFNI_INSTACORR_SEEDBLUR") ;
280      int grad  = (gblur) ? 1.2345f*sblur : 1.0001f*sblur ;
281      MCW_cluster *smask ;
282      float wtsum=1.0f , fac , *qar ;
283      float *sar=(float *)malloc(sizeof(float)*iset->mv->nvals)  ;
284      int qi,qj,qk , ii,ij,ik , qjk,qq , nx,ny,nz,nxy ; register int tt ;
285 
286      if( grad > 0.0f )
287        smask = MCW_spheremask( iset->mv->dx, iset->mv->dy, iset->mv->dz, grad ) ;
288      else
289        smask = MCW_spheremask( 1.0f, 1.0f, 1.0f, -grad ) ;
290 
291      nx = iset->mv->nx ; ny = iset->mv->ny ; nz = iset->mv->nz ; nxy = nx*ny ;
292      ii = ijk % nx ; ik = ijk / nxy ; ij = (ijk-ik*nxy) / nx ;
293 
294      fac = FWHM_TO_SIGMA(sblur); fac = 1.0 / (2.0f * fac*fac); /* for gblur */
295 
296      memcpy(sar,tsar,sizeof(float)*iset->mv->nvals) ;  /* copy of seed */
297 
298      for( kk=1 ; kk < smask->num_pt ; kk++ ){  /* add stuff to this copy */
299        qi  = ii + smask->i[kk] ; if( qi < 0 || qi >= nx ) continue ;
300        qj  = ij + smask->j[kk] ; if( qj < 0 || qj >= ny ) continue ;
301        qk  = ik + smask->k[kk] ; if( qk < 0 || qk >= nz ) continue ;
302        qjk = qi + qj*nx + qk*nxy ;
303        qq  = THD_vectim_ifind( qjk , iset->mv ) ;
304        if( qq >= 0 ){
305          register float wt ;
306          if( gblur ){ float rad=smask->mag[kk]; wt = exp(-fac*rad*rad); } else wt = 1.0f;
307           wtsum += wt ; qar = VECTIM_PTR(iset->mv,qq) ;
308           for( tt=0 ; tt < iset->mv->nvals  ; tt++ ) sar[tt] += wt * qar[tt] ;
309          }
310        }
311        if( wtsum > 1.0f ){  /* usually true, but not always */
312          fac = 1.0f / wtsum ;
313          for( tt=0 ; tt < iset->mv->nvals ; tt++ ) tsar[tt] = fac * sar[tt] ;
314        }
315        free(sar) ; KILL_CLUSTER(smask) ;
316        THD_normalize( iset->mv->nvals , tsar ) ;
317    } /* end of seed blur */
318 
319    RETURN(tsar) ;
320 }
321 
322 /*---------------------------------------------------------------------------*/
323 /*! Compute the instant correlation from spatial voxel index ijk --
324     mri_free() this image when you are done with it!
325 *//*-------------------------------------------------------------------------*/
326 
THD_instacorr(ICOR_setup * iset,int ijk)327 MRI_IMAGE * THD_instacorr( ICOR_setup *iset , int ijk )
328 {
329    int kk ; MRI_IMAGE *qim ; float *qar , *dar , *tsar ; int *ivar ;
330    float sblur=0.0f ;
331    MRI_vectim *mv ;
332    int iter_count , iter , nvals ; float iter_thresh , *nsar=NULL ;
333 
334 ENTRY("THD_instacorr") ;
335 
336 #if 0
337 INFO_message("THD_instacorr: ijk=%d  iset=%p",ijk,(void *)iset) ;
338 if( iset != NULL ) ININFO_message(" iset->mv=%p",(void *)iset->mv) ;
339 #endif
340 
341    if( iset == NULL || iset->mv == NULL || ijk < 0 ) RETURN(NULL) ;
342 
343    /** extract reference time series **/
344 
345    tsar = THD_instacorr_getseed( iset , ijk ) ;
346    if( tsar == NULL ) RETURN(NULL) ; /* bad */
347 
348    /** save seed in iset struct [15 May 2009] **/
349 
350    nvals = iset->mv->nvals ;
351    iset->tseed = (float *)realloc( iset->tseed , sizeof(float)*nvals ) ;
352    memcpy( iset->tseed , tsar , sizeof(float)*nvals ) ;
353 
354    /** do the correlations **/
355 
356    sblur = iset->sblur ;
357 
358    iter_count  = iset->iter_count ; if( iter_count <= 0 ) iter_count = 1 ;
359    iter_thresh = iset->iter_thresh ;
360    if( iter_count > 1 ) nsar = (float *)malloc(sizeof(float)*nvals) ;
361 
362    mv  = (iset->ev == NULL) ? iset->mv : iset->ev ;
363    dar = (float *)malloc(sizeof(float)*iset->mv->nvec) ;
364 
365    for( iter=0 ; iter < iter_count ; iter++ ){  /* iteration: 05 Feb 2015 */
366 
367      switch( iset->cmeth ){  /* the actual work */
368        default:
369        case NBISTAT_PEARSON_CORR:
370          THD_vectim_dotprod ( mv , tsar , dar , 0 ) ; break ;
371 
372        case NBISTAT_SPEARMAN_CORR:
373          THD_vectim_spearman( mv , tsar , dar ) ; break ;
374 
375        case NBISTAT_QUADRANT_CORR:
376          THD_vectim_quadrant( mv , tsar , dar ) ; break ;
377 
378        case NBISTAT_TICTACTOE_CORR:
379          THD_vectim_tictactoe( mv , tsar , dar ) ; break ;
380 
381        case NBISTAT_KENDALL_TAUB:
382          THD_vectim_ktaub( mv , tsar , dar ) ; break ;  /* 29 Apr 2010 */
383 
384        case NBISTAT_BC_PEARSON_M:
385          THD_vectim_pearsonBC( mv,sblur,ijk,0,dar ); break; /* 07 Mar 2011 */
386 
387        case NBISTAT_BC_PEARSON_V:
388          THD_vectim_pearsonBC( mv,sblur,ijk,1,dar ); break; /* 07 Mar 2011 */
389 
390        case NBISTAT_EUCLIDIAN_DIST:/* 04 Apr 2012, ZSS*/
391          THD_vectim_distance( mv , tsar , dar , 0, "inv;n_scale") ; break ;
392 
393        case NBISTAT_CITYBLOCK_DIST:/* 04 Apr 2012, ZSS*/
394          THD_vectim_distance( mv , tsar , dar , 1, "inv;n_scale") ; break ;
395 
396        case NBISTAT_QUANTILE_CORR: /* 11 May 2012 */
397          THD_vectim_quantile( mv , tsar , dar ) ; break ;
398      }
399 
400      if( iter+1 < iter_count ){  /* threshold, compute new seed [05 Feb 2015] */
401        int ii , ngood=0 ; float *vv ;
402        for( ii=0 ; ii < nvals ; ii++ ) nsar[ii] = 0.0f ;
403 
404        for( kk=0 ; kk < mv->nvec ; kk++ ){
405          vv = VECTIM_PTR(mv,kk) ;
406          if( dar[kk] >= iter_thresh ){
407            ngood++ ; for( ii=0 ; ii < nvals ; ii++ ) nsar[ii] += vv[ii] ;
408          } else if( dar[kk] <= -iter_thresh ){
409            ngood++ ; for( ii=0 ; ii < nvals ; ii++ ) nsar[ii] -= vv[ii] ;
410          }
411        }
412        if( ngood == 0 ){  /* nothing above threshold?! */
413          WARNING_message("InstaCorr -- nothing above threshold %g on iteration #%d",
414                          iter_thresh , iter+1 ) ;
415          break ;
416        }
417 #if 0
418 ININFO_message("InstaCorr: iteration #%d has %d supra-threshold",iter+1,ngood) ;
419 #endif
420        THD_normalize( nvals , nsar ) ;
421        memcpy( tsar , nsar , sizeof(float)*nvals ) ;
422 
423      } else if( nsar != NULL ){
424        free(nsar) ; nsar = NULL ;
425      }
426 
427    } /* end of iterations */
428 
429    /** put them into the output image **/
430 
431    qim  = mri_new_vol( mv->nx , mv->ny , mv->nz , MRI_float ) ;
432    qar  = MRI_FLOAT_PTR(qim) ;
433    ivar = mv->ivec ;
434    for( kk=0 ; kk < mv->nvec ; kk++ ) qar[ivar[kk]] = dar[kk] ;
435 
436    /** e finito **/
437 
438    free(dar) ; free(tsar) ; RETURN(qim) ;
439 }
440 
441 /*---------------------------------------------------------------------------*/
442 /*! Compute the instant correlation from tsar as the seed series, from time
443     indexes ibot to itop-1 -- mri_free() result when you are done with it
444 *//*-------------------------------------------------------------------------*/
445 
THD_instacorr_section(ICOR_setup * iset,float * tsar,int ibot,int itop)446 MRI_IMAGE * THD_instacorr_section( ICOR_setup *iset, float *tsar, int ibot, int itop )
447 {
448    int kk ; MRI_IMAGE *qim ; float *qar , *dar ; int *ivar ;
449    MRI_vectim *mv ;
450 
451 ENTRY("THD_instacorr_section") ;
452 
453    if( iset == NULL || iset->mv == NULL || tsar == NULL ) RETURN(NULL) ;
454    if( ibot < 0 ) ibot = 0 ;
455    if( itop > iset->mv->nvals-1 ) itop = iset->mv->nvals-1 ;
456    if( itop-ibot < 5 ) RETURN(NULL) ;  /* too short to correlate */
457 
458 #if 0
459 ININFO_message("section: ibot=%d itop=%d",ibot,itop) ;
460 #endif
461 
462    /** do the correlations **/
463 
464    dar = (float *)malloc(sizeof(float)*iset->mv->nvec) ;
465 
466    mv = (iset->ev == NULL) ? iset->mv : iset->ev ;
467 
468    switch( iset->cmeth ){
469      default:
470      case NBISTAT_PEARSON_CORR:
471        THD_vectim_pearson_section( mv , tsar , dar , ibot,itop ) ; break ;
472 
473 #if 0
474      case NBISTAT_SPEARMAN_CORR:
475        THD_vectim_spearman( mv , tsar , dar ) ; break ;
476 
477      case NBISTAT_QUADRANT_CORR:
478        THD_vectim_quadrant( mv , tsar , dar ) ; break ;
479 
480      case NBISTAT_TICTACTOE_CORR:
481        THD_vectim_tictactoe( mv , tsar , dar ) ; break ;
482 
483      case NBISTAT_KENDALL_TAUB:
484        THD_vectim_ktaub( mv , tsar , dar ) ; break ;  /* 29 Apr 2010 */
485 
486      case NBISTAT_BC_PEARSON_M:
487        THD_vectim_pearsonBC( mv,sblur,ijk,0,dar ); break; /* 07 Mar 2011 */
488 
489      case NBISTAT_BC_PEARSON_V:
490        THD_vectim_pearsonBC( mv,sblur,ijk,1,dar ); break; /* 07 Mar 2011 */
491 
492      case NBISTAT_EUCLIDIAN_DIST:/* 04 Apr 2012, ZSS*/
493        THD_vectim_distance( mv , tsar , dar , 0, "inv;n_scale") ; break ;
494 
495      case NBISTAT_CITYBLOCK_DIST:/* 04 Apr 2012, ZSS*/
496        THD_vectim_distance( mv , tsar , dar , 1, "inv;n_scale") ; break ;
497 
498      case NBISTAT_QUANTILE_CORR: /* 11 May 2012 */
499        THD_vectim_quantile( mv , tsar , dar ) ; break ;
500 #endif
501    }
502 
503    /** put them into the output image **/
504 
505    qim  = mri_new_vol( mv->nx , mv->ny , mv->nz , MRI_float ) ;
506    qar  = MRI_FLOAT_PTR(qim) ;
507    ivar = mv->ivec ;
508    for( kk=0 ; kk < mv->nvec ; kk++ ) qar[ivar[kk]] = dar[kk] ;
509 
510    /** e finito **/
511 
512    free(dar) ; RETURN(qim) ;
513 }
514 
515 /*---------------------------------------------------------------------------*/
516 
THD_instacorr_collection(ICOR_setup * iset,int ijk)517 MRI_IMARR * THD_instacorr_collection( ICOR_setup *iset , int ijk )
518 {
519    int ibot , itop , nii ,  kk ;
520    float *tsar ;
521    MRI_IMARR *outar ; MRI_IMAGE *outim ;
522 
523    /** extract reference time series **/
524 
525    tsar = THD_instacorr_getseed( iset , ijk ) ;
526    if( tsar == NULL ) RETURN(NULL) ; /* bad */
527 
528    INIT_IMARR(outar) ;
529 
530    for( nii=0,ibot=0 ; nii < iset->cnum ; ibot+=iset->cstep,nii++ ){
531      itop = ibot + iset->clen-1 ;
532      outim = THD_instacorr_section( iset , tsar , ibot,itop ) ;
533      ADDTO_IMARR(outar,outim) ;
534    }
535    free(tsar) ;
536    RETURN(outar) ;
537 }
538