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