1 static int verb = 1 ;
2
3 static void NCO_help(void) ; /* prototype */
4
5 /*---------------------------------------------------------------------------*/
6 #ifndef FLOATIZE
7 # include "matrix.h" /* double precision */
8 # define MTYPE double
9 # define NI_MTYPE NI_DOUBLE
10 # define QEPS 1.e-6
11 #else
12 # include "matrix_f.h" /* single precision */
13 # define MTYPE float
14 # define NI_MTYPE NI_FLOAT
15 # define QEPS 1.e-4
16 #endif
17
18 #include "mrilib.h" /* Keep after decision about matrix.h inclusion
19 ZSS Nov. 21 2014*/
20
21 /*---------------------------------------------------------------------------*/
22 /* Everything needed to solve one linear problem (baseline and full models) */
23
24 typedef struct {
25 matrix x_all ; /* the whole matrix (before censoring): nall X p */
26 matrix x_full ; /* after censoring: nfull X p */
27
28 matrix xtxinv_full ; /* p X p */
29 matrix xtxinvxt_full ; /* p X nfull */
30
31 matrix x_base ; /* nfull X q */
32 matrix xtxinvxt_base ; /* q X nfull */
33 } linear_setup ;
34
35 #define INIT_linear_setup(ls) \
36 do{ matrix_initialize(&((ls).x_all)) ; \
37 matrix_initialize(&((ls).x_full)) ; \
38 matrix_initialize(&((ls).xtxinv_full)) ; \
39 matrix_initialize(&((ls).xtxinvxt_full)) ; \
40 matrix_initialize(&((ls).x_base) ; \
41 matrix_initialize(&((ls).xtxinvxt_base) ; \
42 } while(0)
43
44 /*---------------------------------------------------------------------------*/
45 /* Specify one HRF model function */
46
47 typedef struct {
48 int parnum ; /* number of parameters */
49 float tbot , ttop ; /* start and stop times */
50 float *parval , *parbot , *partop ; /* parameter values and ranges */
51 void *qpar ; /* other information */
52 float (*f)(float,int,float *,void *); /* function to evaluate HRF */
53 } hrf_model ;
54
55 /*---------------------------------------------------------------------------*/
56 /* Everything needed to solve the nonlinear problem. */
57
58 typedef struct {
59 int polort, /* number of polynomial parameters PER run */
60 qp, /* total number of polynomial baseline parameters */
61 q , /* total number of baseline parameters */
62 p , /* total number of linear parameters (matrix columns) */
63 nall , /* number of input data rows (before censoring) */
64 nfull ; /* number of input data rows (after censoring) */
65
66 int ngood ; /* 0 if no censoring, or nfull if censoring */
67 int *goodlist ; /* goodlist[i] = index of i-th good row in all data */
68 int run_num ; /* number of runs in total data */
69 int *run_start; /* run_start[i] = data index of start of run #i */
70 float tr ; /* time step of data */
71
72 int nprob ; /* number of linear setups (slices) */
73 linear_setup *prob ; /* one per slice */
74 float *prob_toff ; /* time shift for each prob */
75
76 /* Data time series from input dataset,
77 separated into slices and in time-then-space order: */
78 int *nvec ; /* nvec[i] = # vectors (length=nfull) in prob #i */
79 float **vec ; /* vec[i] = ptr to nvec[i] data vectors */
80 int **ivec ; /* ivec[i][k] = 3D index of where k-th data vector is from */
81 int nx,ny,nz; /* spatial dimensions of input grid */
82
83 int num_hrf ; /* number of HRF model functions we're finding */
84 hrf_model *hrf ; /* description of HRF model function */
85
86 int num_times ; /* number of -stimtime options */
87 char **times_label ; /* label for each option */
88 MRI_IMAGE **times ; /* times for each option */
89 int *times_hrfind ; /* which HRF model function to use */
90 } nonlinear_setup ;
91
92 /*---------------------------------------------------------------------------*/
93 /* Given the x_all matrix, finish setting up the rest of the stuff
94 (i.e., extract the full and base matrices, compute pseudo-inverses).
95 -----------------------------------------------------------------------------*/
96
NCO_finish_linear_setup(linear_setup * ls,int ngd,int * gdlist,int qb)97 void NCO_finish_linear_setup( linear_setup *ls, int ngd, int *gdlist, int qb )
98 {
99 matrix *xf ; /* stand-in for full matrix */
100
101 ENTRY("NCO_finish_linear_setup") ;
102
103 if( !ISVALID_MATRIX(ls->x_all) ){
104 WARNING_message("Unprepared call to NCO_finish_linear_setup"); EXRETURN;
105 }
106
107 /* on each entry, x_all was changed in the HRF part,
108 so need to re-extract and pseudo-invert the x_full sub-matrix */
109
110 if( ISVALID_MATRIX(ls->x_full) ) matrix_destroy( &(ls->x_full) ) ;
111
112 if( ngd > 0 && gdlist != NULL ){ /* if any censoring was done */
113 xf = &(ls->x_full) ;
114 matrix_extract_rows( ls->x_all , ngd,gdlist , xf ) ;
115 } else { /* no censoring was done */
116 xf = &(ls->x_all) ;
117 }
118 matrix_psinv( *xf , &(ls->xtxinv_full) , &(ls->xtxinvxt_full) ) ;
119
120 /* only on first entry should baseline part of x_all change */
121
122 if( qb > 0 && !ISVALID_MATRIX(ls->x_base) ){
123 int *qlist=malloc(sizeof(int)*qb) , ii ;
124 for( ii=0 ; ii < qb ; ii++ ) qlist[ii] = ii ;
125 matrix_extract_cols( *xf , qb,qlist , &(ls->x_base) ) ;
126 free((void *)qlist) ;
127 matrix_psinv( ls->x_base , NULL , &(ls->xtxinvxt_base) ) ;
128 }
129
130 EXRETURN ;
131 }
132
133 /*---------------------------------------------------------------------------*/
134 /* Erase the forgettable parts of a linear_setup struct. */
135
NCO_clear_linear_setup(linear_setup * ls,int dobase)136 void NCO_clear_linear_setup( linear_setup *ls , int dobase )
137 {
138 matrix_destroy( &(ls->x_full) ) ;
139 matrix_destroy( &(ls->xtxinv_full) ) ;
140 matrix_destroy( &(ls->xtxinvxt_full) ) ;
141 if( dobase ){
142 matrix_destroy( &(ls->x_base) ) ;
143 matrix_destroy( &(ls->xtxinvxt_base) ) ;
144 }
145 return ;
146 }
147
148 /*---------------------------------------------------------------------------*/
149 /* Extract all slice data in the mask from the given dataset.
150 Output is an image of time series vectors, and an image of their
151 3D indexes in the dataset array.
152 -----------------------------------------------------------------------------*/
153
NCO_extract_slice(THD_3dim_dataset * dset,int ss,byte * mask,int ngood,int * goodlist)154 MRI_IMARR * NCO_extract_slice( THD_3dim_dataset *dset, int ss,
155 byte *mask, int ngood, int *goodlist )
156 {
157 int ii,jj,kk,vv,uu,ww , nx,ny,nz , nt,ng , nin;
158 MRI_IMAGE *sim ; float *sar , *tar , *qar ;
159 MRI_IMAGE *iim ; int *iar ;
160 MRI_IMARR *imar ;
161
162 ENTRY("NCO_extract_slice") ;
163
164 nx = DSET_NX(dset); ny = DSET_NY(dset);
165 nz = DSET_NZ(dset); nt = DSET_NVALS(dset);
166 if( ss < 0 || ss >= nz ) ERROR_exit("illegal slice index %d!?",ss) ;
167
168 /* number of output time points */
169 ng = (ngood > 0 && goodlist != NULL) ? ngood : nt ;
170
171 vv = ss*nx*ny ; /* slice index offset into 3D arrays */
172 if( mask != NULL ){ /* count number of voxels in mask & in this slice */
173 for( nin=jj=0 ; jj < ny ; jj++ ){
174 uu = vv + jj*nx ;
175 for( ii=0 ; ii < nx ; ii++ ) if( mask[ii+uu] ) nin++ ;
176 }
177 if( nin == 0 ) RETURN(NULL) ; /* nothing in the mask */
178 } else {
179 nin = nx*ny ; /* use all voxels in slice */
180 }
181
182 /* create output image */
183
184 sim = mri_new(ng,nin,MRI_float) ; sar = MRI_FLOAT_PTR(sim) ;
185 iim = mri_new(nin,1 ,MRI_int) ; iar = MRI_INT_PTR (iim) ;
186
187 /* perhaps need temp staging area for array from dataset? */
188
189 if( ng < nt ) tar = (float *)malloc(sizeof(float)*nt) ;
190 else tar = NULL ;
191
192 /* loop over voxels in this slice */
193
194 for( kk=jj=0 ; jj < ny ; jj++ ){ /* kk=output voxel index */
195 uu = vv + jj*nx ; /* uu=index to jj-th row in slice */
196 for( ii=0 ; ii < nx ; ii++ ){
197 if( mask != NULL && mask[ii+uu] == 0 ) continue ; /* skip this'n */
198 qar = sar + kk*ng ; /* where to store output */
199 if( tar != NULL ){ /* get temp array, then sub-sample */
200 THD_extract_array( ii+uu , dset , 0 , tar ) ;
201 for( ww=0 ; ww < ng ; ww++ ) qar[ww] = tar[goodlist[ww]] ;
202 } else { /* get directly into output */
203 THD_extract_array( ii+uu , dset , 0 , qar ) ;
204 }
205 iar[kk] = ii+uu ; /* save index of where this data vector came from */
206 kk++ ; /* increment count of saved voxel vectors */
207 }
208 }
209
210 if( tar != NULL ) free((void *)tar) ; /* toss some trash */
211
212 INIT_IMARR(imar); ADDTO_IMARR(imar,sim); ADDTO_IMARR(imar,iim); RETURN(imar);
213 }
214
215 /*===========================================================================*/
216 /*---------------------- Main program (not much here yet) -------------------*/
217
main(int argc,char * argv[])218 int main( int argc , char *argv[] )
219 {
220 int iarg=1 ;
221
222 THD_3dim_dataset *inset=NULL , *maskset=NULL ;
223 int automask=0 ;
224 float tr=0.0f ;
225 nonlinear_setup *nls ;
226 MRI_IMARR *baseim ;
227 MRI_IMAGE *concatim=NULL ;
228 int num_CENSOR=0 ;
229 int_triple *abc_CENSOR=NULL ;
230 char *tpat=NULL ;
231
232 int ii,jj,kk , nslice , nt , nvectot ;
233 int nmask=0 ; byte *mask=NULL ;
234 MRI_IMAGE *sim,*iim ; MRI_IMARR *imar ; float *sar ; int *iar ;
235
236 /*----- help the pitiful user? -----*/
237
238 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ NCO_help(); exit(0); }
239
240 #ifdef USING_MCW_MALLOC
241 enable_mcw_malloc() ;
242 #endif
243 PRINT_VERSION("3dNeocon") ; AUTHOR("RW Cox");
244 mainENTRY("3dNeocon main") ; machdep() ;
245 AFNI_logger("3dNeocon",argc,argv) ;
246
247 /*----- read arguments -----*/
248
249 nls = (nonlinear_setup *)calloc(1,sizeof(nonlinear_setup)) ;
250 nls->polort = -666 ; /* automatic */
251 INIT_IMARR(baseim) ;
252
253 while( iarg < argc ){
254
255 /*----------*/
256
257 if( strcmp(argv[iarg],"-input") == 0 ){
258 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]) ;
259 if( inset != NULL )
260 ERROR_exit("Can't have two -input arguments!") ;
261 inset = THD_open_dataset( argv[++iarg] ) ;
262 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
263 if( DSET_NVALS(inset) < 9 )
264 ERROR_exit("Input dataset has less than 9 time points!") ;
265 iarg++ ; continue ;
266 }
267
268 /*----------*/
269
270 if( strcmp(argv[iarg],"-mask") == 0 ){
271 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]) ;
272 if( maskset != NULL )
273 ERROR_exit("Can't have two -mask arguments!") ;
274 if( automask )
275 ERROR_exit("Can't combine -mask and -automask!") ;
276 maskset = THD_open_dataset( argv[++iarg] ) ;
277 CHECK_OPEN_ERROR(maskset,argv[iarg]) ;
278 iarg++ ; continue ;
279 }
280
281 /*----------*/
282
283 if( strcmp(argv[iarg],"-automask") == 0 ){
284 if( maskset != NULL ) ERROR_exit("Can't combine -automask and -mask!");
285 automask = 1 ;
286 iarg++ ; continue ;
287 }
288
289 /*----------*/
290
291 if( strcmp(argv[iarg],"-TR") == 0 ){
292 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]) ;
293 tr = (float)strtod(argv[++iarg],NULL) ;
294 if( tr <= 0.0f ) ERROR_exit("Can't read valid value after -TR") ;
295 iarg++ ; continue ;
296 }
297
298 /*----------*/
299
300 if( strcmp(argv[iarg],"-polort") == 0 ){
301 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]) ;
302 iarg++ ;
303 if( argv[iarg][0] == 'A' ) nls->polort = -1 ;
304 else nls->polort = (int)strtod(argv[iarg],NULL) ;
305 iarg++ ; continue ;
306 }
307
308 /*----------*/
309
310 if( strcmp(argv[iarg],"-baseline") == 0 ){
311 MRI_IMAGE *bim ;
312 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]) ;
313 bim = mri_read_1D( argv[++iarg] ) ;
314 if( bim == NULL ) ERROR_exit("Can't read -baseline '%s'",argv[iarg]) ;
315 ADDTO_IMARR(baseim,bim) ;
316 iarg++ ; continue ;
317 }
318
319 /*----------*/
320
321 if( strcmp(argv[iarg],"-concat") == 0 ){
322 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]);
323 if( concatim != NULL ) ERROR_exit("Can't have 2 -concat arguments!");
324 concatim = mri_read_1D( argv[++iarg] ) ;
325 if( concatim == NULL ) ERROR_exit("Can't read -concat '%s'",argv[iarg]);
326 iarg++ ; continue ;
327 }
328
329 /*----------*/
330
331 if( strncmp(argv[iarg],"-CENSOR",7) == 0 || /* adapted from */
332 strncmp(argv[iarg],"-censorTR",9) == 0 ){ /* 3dDeconvolve */
333
334 NI_str_array *nsar ;
335 char *src=malloc(1), *cpt, *dpt ;
336 int ns, r,a,b, nerr=0 ; int_triple rab ;
337
338 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]);
339
340 *src = '\0' ; /* cat all following options until starts with '-' */
341 for( iarg++ ; iarg < argc && argv[iarg][0] != '-' ; iarg++ ){
342 ns = strlen(argv[iarg]) ; if( ns == 0 ) continue ;
343 src = realloc(src,strlen(src)+ns+2) ;
344 strcat(src," ") ; strcat(src,argv[iarg]) ;
345 }
346 if( *src == '\0' ) ERROR_exit("Bad argument after -CENSORTR") ;
347 nsar = NI_decode_string_list( src , "," ) ; /* break into substrings */
348 for( ns=0 ; ns < nsar->num ; ns++ ){ /* loop over substrings */
349 cpt = nsar->str[ns] ; dpt = strchr(cpt,':') ; r = 0 ;
350 if( *cpt == '\0' ) continue ; /* skip an empty string */
351 if( dpt != NULL ){ /* found 'run:' */
352 if( *cpt == '*' ){ /* wildcard = all runs */
353 r = -666 ;
354 } else {
355 r = (int)strtol(cpt,NULL,10) ;
356 if( r <= 0 ){ /* skip out */
357 ERROR_message("-CENSORTR %s -- run index '%d' is bad!",nsar->str[ns],r);
358 nerr++ ; continue ;
359 }
360 }
361 cpt = dpt+1 ; /* skip to character after ':' */
362 if( *cpt == '\0' ){ /* skip out */
363 ERROR_message("-CENSORTR %s -- no data after run index!",nsar->str[ns]);
364 nerr++ ; continue ;
365 }
366 }
367 a = (int)strtol(cpt,&dpt,10) ; /* get first index number */
368 if( a < 0 ){ /* skip out */
369 ERROR_message("-CENSORTR %s -- time index '%d' is bad!",nsar->str[ns],a);
370 nerr++ ; continue ;
371 }
372 if( *dpt == '\0' ){ /* no second number */
373 b = a ;
374 } else { /* get second number */
375 for( dpt++ ; *dpt != '\0' && !isdigit(*dpt) ; dpt++ ) ; /*nada*/
376 b = (int)strtol(dpt,NULL,10) ;
377 if( b < a || b < 0 ){ /* skip out */
378 ERROR_message("-CENSORTR %s -- time indexes '%d' to '%d' is bad!",
379 nsar->str[ns],a,b);
380 nerr++ ; continue ;
381 }
382 }
383 abc_CENSOR = (int_triple *)realloc( abc_CENSOR ,
384 sizeof(int_triple)*(num_CENSOR+1) );
385 rab.i = r; rab.j = a; rab.k = b; abc_CENSOR[num_CENSOR++] = rab ;
386 } /* end of loop over -CENSORTR strings */
387 if( nerr > 0 ) ERROR_exit("Can't proceed after -CENSORTR errors!") ;
388 NI_delete_str_array(nsar) ; free(src) ;
389 continue ; /* next option */
390 }
391
392 /*----------*/
393
394 if( strcasecmp(argv[iarg],"-verb") == 0 ){ verb++ ; iarg++ ; continue ; }
395 if( strcasecmp(argv[iarg],"-quiet")== 0 ){ verb=0 ; iarg++ ; continue ; }
396
397 /*----------*/
398
399 if( strcmp(argv[iarg],"-tpattern") == 0 ){
400 if( iarg == argc-1 ) ERROR_exit("Need argument after '%s'",argv[iarg]);
401 tpat = argv[++iarg] ;
402 iarg++ ; continue ;
403 }
404
405 /*----------*/
406
407 if( strcmp(argv[iarg],"-stimtime") == 0 ){ iarg++ ; continue ; }
408 if( strcmp(argv[iarg],"-threshtype") == 0 ){ iarg++ ; continue ; }
409 if( strcmp(argv[iarg],"-fitts") == 0 ){ iarg++ ; continue ; }
410 if( strcmp(argv[iarg],"-errts") == 0 ){ iarg++ ; continue ; }
411 if( strcmp(argv[iarg],"-cbucket") == 0 ){ iarg++ ; continue ; }
412
413 /*==========*/
414
415 ERROR_exit("Unknown argument: '%s'",argv[iarg]) ;
416
417 } /*----- end of scanning args -----*/
418
419 /*--- check for input errors or inconsistencies ---*/
420
421 if( iarg < argc ) WARNING_message("arguments after last option???") ;
422
423 if( inset == NULL ) ERROR_exit("No -input dataset given?") ;
424
425 /* check CENSOR command for run indexes: should all have them, or none */
426
427 if( abc_CENSOR != NULL ){
428 int ic , rr , nzr=0 ;
429 for( ic=0 ; ic < num_CENSOR ; ic++ ) /* count number with run != 0 */
430 if( abc_CENSOR[ic].i ) nzr++ ;
431 if( nzr > 0 && nzr < num_CENSOR )
432 WARNING_message("%d -CENSORTR commands have run: numbers and %d do not!" ,
433 nzr , num_CENSOR-nzr ) ;
434 }
435
436 /*--- setup timing information ---*/
437
438 if( tr <= 0.0f ){
439 tr = DSET_TR(inset) ;
440 if( tr <= 0.0f ){
441 WARNING_message("no TR in dataset and no -TR option ==> TR = 1 s") ;
442 tr = 1.0f ;
443 }
444 } else {
445 float dd = DSET_TR(inset) ;
446 if( dd > 0.0f && fabsf(dd-tr) > 0.001f )
447 WARNING_message("-TR %g overrides dataset TR=%g",tr,dd) ;
448 }
449 nls->tr = tr ;
450
451 if( tpat != NULL ){
452 nls->prob_toff = TS_parse_tpattern( DSET_NZ(inset), nls->tr , tpat ) ;
453 if( nls->prob_toff == NULL )
454 nls->prob_toff = (float *)calloc(sizeof(float),DSET_NZ(inset)) ;
455 } else {
456 nls->prob_toff = (float *)calloc(sizeof(float),DSET_NZ(inset)) ;
457 if( inset->taxis != NULL && inset->taxis->toff_sl != NULL )
458 memcpy( nls->prob_toff, inset->taxis->toff_sl,
459 sizeof(float)*DSET_NZ(inset) ) ;
460 }
461
462 /*----- setup run start indexes -----*/
463
464 nt = DSET_NVALS(inset) ;
465
466 if( DSET_IS_TCAT(inset) ){ /* from dataset catenated on command line */
467
468 if( concatim != NULL )
469 WARNING_message("-concat ignored since -input catenated on command line");
470 nls->run_num = inset->tcat_num ;
471 nls->run_start = (int *)malloc( sizeof(int)*nls->run_num ) ;
472 nls->run_start[0] = 0 ;
473 for( ii=0 ; ii < nls->run_num-1 ; ii++ )
474 nls->run_start[ii+1] = nls->run_start[ii] + inset->tcat_len[ii] ;
475
476 } else if( concatim != NULL ){ /* from -concat */
477
478 float *car = MRI_FLOAT_PTR(concatim) ;
479 nls->run_num = concatim->nvox ;
480 nls->run_start = (int *)malloc( sizeof(int)*nls->run_num ) ;
481 for( ii=0 ; ii < nls->run_num ; ii++ ) nls->run_start[ii] = (int)car[ii] ;
482 mri_free(concatim); concatim = NULL;
483 jj = 0 ;
484 if( nls->run_start[0] != 0 ) jj++ ;
485 for( ii=1 ; ii < nls->run_num ; ii++ )
486 if( nls->run_start[ii] <= nls->run_start[ii-1] ) jj++ ;
487 if( jj > 0 ) ERROR_exit("Disordered run start indexes in -concat") ;
488 if( nls->run_start[nls->run_num-1] >= nt )
489 ERROR_exit("last -concat value %d is after end of input dataset (%d)?!",
490 nls->run_start[nls->run_num-1] , nt-1 ) ;
491
492 } else { /* default = just 1 big happy run */
493
494 nls->run_num = 1 ;
495 nls->run_start = (int *)malloc(sizeof(int)) ;
496 nls->run_start[0] = 0 ;
497 }
498
499 /*--- estimate desirable polort from max block duration ---*/
500
501 { int ibot, itop, ilen, lmax=0 ; float dmax ;
502 for( ii=0 ; ii < nls->run_num ; ii++ ){
503 ibot = nls->run_start[ii] ;
504 itop = (ii < nls->run_num-1) ? nls->run_start[ii+1] : nt ;
505 ilen = itop - ibot ; lmax = MAX(lmax,ilen) ;
506 }
507 dmax = tr * lmax ; ilen = 1 + (int)floor(dmax/150.0) ;
508 if( nls->polort < -1 ){
509 nls->polort = ilen ;
510 INFO_message("longest run=%.1f s; automatic polort=%d",dmax,ilen) ;
511 } else if( nls->polort >= 0 && ilen > nls->polort ){
512 WARNING_message(
513 "input polort=%d; longest run=%.1fs; recommended polort=%d",
514 nls->polort , dmax , ilen ) ;
515 }
516 }
517
518 /*----- mask-ification -----*/
519
520 if( automask ){
521 if( verb > 1 ) INFO_message("Making automask from input dataset") ;
522 mask = THD_automask( inset ) ;
523 if( mask == NULL ) ERROR_exit("automask-ing fails?!") ;
524 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
525 if( verb ) INFO_message("%d voxels in automask",nmask) ;
526 } else if( maskset != NULL ){
527 if( DSET_NVOX(maskset) != DSET_NVOX(inset) )
528 ERROR_exit("-mask and -input datasets not the same size!") ;
529 if( !EQUIV_GRIDS(maskset,inset) )
530 WARNING_message("-mask and -input datasets have differing grids!") ;
531 if( verb > 1 ) INFO_message("reading mask dataset") ;
532 mask = THD_makemask( maskset , 0 , 1.0f,-1.0f ) ;
533 DSET_delete(maskset) ;
534 if( mask == NULL ) ERROR_exit("mask is empty?!") ;
535 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
536 if( verb ) INFO_message("%d voxels in mask",nmask) ;
537 } else {
538 nmask = DSET_NVOX(inset) ; mask = NULL ;
539 if( verb ) INFO_message("no masking ==> will process all %d voxels",nmask) ;
540 }
541
542 /*----- construct the goodlist (from the censoring data) -----*/
543
544 nls->nall = nt ;
545 if( num_CENSOR == 0 || abc_CENSOR == NULL ){
546 nls->ngood = 0 ; nls->goodlist = NULL ; nls->nfull = nt ;
547 if( verb > 1 )
548 INFO_message("no censoring ==> all %d time points used",nt) ;
549 } else {
550 int ic,it , rr , aa,bb , nerr=0 , bbot,btop , nblk=nls->run_num ;
551 byte *cenar=(byte *)malloc(sizeof(byte)*nt) ;
552 for( it=0 ; it < nt ; it++ ) cenar[it] = 1 ;
553 for( ic=0 ; ic < num_CENSOR ; ic++ ){ /* loop over CENSOR commands */
554 rr = abc_CENSOR[ic].i ;
555 aa = abc_CENSOR[ic].j ; if( aa < 0 ) continue ; /* shouldn't happen */
556 bb = abc_CENSOR[ic].k ; if( bb < aa ) continue ; /* shouldn't happen */
557 if( rr == -666 ){ /* run = wildcard ==> expand to nblk new triples */
558 int_triple rab ;
559 abc_CENSOR = (int_triple *)realloc( abc_CENSOR ,
560 sizeof(int_triple)*(num_CENSOR+nblk) );
561 for( rr=1 ; rr <= nblk ; rr++ ){
562 rab.i = rr; rab.j = aa; rab.k = bb; abc_CENSOR[num_CENSOR++] = rab;
563 }
564 continue ; /* skip to next one */
565 }
566 if( rr > 0 ){ /* convert local indexes to global */
567 if( rr > nblk ){ /* stupid user */
568 ERROR_message("-CENSORTR %d:%d-%d has run index out of range 1..%d",
569 rr,aa,bb , nblk ) ;
570 nerr++ ; aa = -66666666 ;
571 } else {
572 bbot = nls->run_start[rr-1] ; /* start index of block #rr */
573 btop = (rr < nblk) ? nls->run_start[rr]-1 : nt-1 ; /* last index */
574 if( aa+bbot > btop ){ /* WTF? */
575 WARNING_message(
576 "-CENSORTR %d:%d-%d has start index past end of run (%d) - IGNORING",
577 rr,aa,bb,btop-bbot ) ; aa = -66666666 ;
578 } else if( bb+bbot > btop ){ /* oopsie */
579 WARNING_message(
580 "-CENSORTR %d:%d-%d has stop index past end of run (%d) - STOPPING THERE",
581 rr,aa,bb,btop-bbot ) ;
582 }
583 aa += bbot ; bb += bbot ; if( bb > btop ) bb = btop ;
584 }
585 } else { /* global indexes: check for stupidities */
586 if( aa >= nt ){
587 WARNING_message(
588 "-CENSORTR %d..%d has start index past end of data (%d) - IGNORING",
589 rr,aa,bb,nt-1 ) ; aa = -66666666 ;
590 } else if( bb > nt ){
591 WARNING_message(
592 "-CENSORTR %d..%d has stop index past end of data (%d) - STOPPING THERE",
593 rr,aa,bb,nt-1 ) ; bb = nt-1 ;
594 }
595 }
596 if( aa < 0 || aa >= nt ) continue ; /* nothing to do */
597 if( bb < aa || bb >= nt ) continue ;
598 if( verb > 1 )
599 ININFO_message("applying -CENSORTR global time indexes %d..%d",aa,bb) ;
600 for( it=aa ; it <= bb ; it++ ) cenar[it] = 0 ;
601 } /* end of loop over CENSOR commands */
602 if( nerr > 0 ) ERROR_exit("Can't continue! Fix the -CENSORTR error%s",
603 (nerr==1) ? "." : "s." ) ;
604 free((void *)abc_CENSOR) ; abc_CENSOR = NULL ; num_CENSOR = 0 ;
605 nls->goodlist = (int *)malloc(sizeof(int)*nt) ;
606 for( ic=it=0 ; it < nt ; it++ ){
607 if( cenar[it] ) nls->goodlist[ic++] = it ;
608 }
609 nls->ngood = nls->nfull = ic ;
610 if( ic < 3 ) ERROR_exit("Number of points surviving censoring = %d",ic) ;
611 free((void *)cenar) ;
612 if( verb > 1 )
613 INFO_message("censoring ==> %d time points used (out of %d total)",
614 nls->ngood , nt ) ;
615 }
616
617 /*----- extract time series from input dataset -----*/
618
619 if( verb > 1 ) INFO_message("loading time series vectors") ;
620
621 DSET_load(inset); CHECK_LOAD_ERROR(inset);
622
623 nls->nz = DSET_NZ(inset); nls->ny = DSET_NY(inset); nls->nx = DSET_NX(inset);
624 nls->nprob = nls->nz ;
625 nls->nvec = (int *) calloc(nls->nprob,sizeof(int) ) ;
626 nls->vec = (float **)calloc(nls->nprob,sizeof(float *)) ;
627 nls->ivec = (int **) calloc(nls->nprob,sizeof(int *) ) ;
628 for( nvectot=kk=0 ; kk < nls->nz ; kk++ ){
629 imar = NCO_extract_slice( inset, kk, mask, nls->ngood,nls->goodlist ) ;
630 if( imar == NULL || IMARR_COUNT(imar) != 2 ){
631 WARNING_message("No data extracted in slice #%d",kk) ;
632 nls->vec[kk] = NULL ; nls->ivec[kk] = NULL ; nls->nvec[kk] = 0 ;
633 } else {
634 sim = IMARR_SUBIM(imar,0) ; sar = MRI_FLOAT_PTR(sim) ;
635 iim = IMARR_SUBIM(imar,1) ; iar = MRI_INT_PTR (iim) ;
636 nls->vec[kk] = MRI_FLOAT_PTR(sim) ; mri_clear_data_pointer(sim) ;
637 nls->ivec[kk] = MRI_INT_PTR(iim) ; mri_clear_data_pointer(iim) ;
638 nls->nvec[kk] = iim->nx ; nvectot += iim->nx ;
639 DESTROY_IMARR(imar) ;
640 }
641 }
642 DSET_unload(inset) ;
643 if( nvectot == 0 ) ERROR_exit("no time series to analyze?!") ;
644 else if( verb ) INFO_message("analyzing %d time series",nvectot) ;
645
646 /************/
647 exit(0) ;
648 }
649
650 /*---------------------------------------------------------------------------*/
651
NCO_help(void)652 static void NCO_help(void)
653 {
654 printf("\n"
655 "This program fits a mixed linear/nonlinear deconvolution model to FMRI\n"
656 "time series. At present, it is experimental, so be careful out there.\n"
657 "\n"
658 "Usage: 3dNeocon ...\n"
659 "\n"
660 "Data Arguments\n"
661 "--------------\n"
662 "These arguments specify the input data and things about it.\n"
663 "\n"
664 " -input dname = read 3D+time dataset 'dname'\n"
665 " -mask mname = read dataset 'mname' as a mask\n"
666 " -automask = compute a mask from the 3D+time dataset\n"
667 " -TR tt = use 'tt' as the TR (in seconds)\n"
668 " -CENSORTR clist = like in 3dDeconvolve\n"
669 " -concat rname = like in 3dDeconvolve\n"
670 " -tpattern ppp = set the slice timing pattern to 'ppp', as in\n"
671 " 3dTshift, to override whatever slice timing\n"
672 " information is in the '-input' dataset header;\n"
673 " in particular, use 'zero' or 'equal' to specify\n"
674 " that all slices are to be treated as acquired\n"
675 " simultaneously (e.g., 3D imaging or 2D imaging\n"
676 " with slice timing correction already performed)\n"
677 "\n"
678 "Baseline Model Arguments\n"
679 "------------------------\n"
680 "These arguments set up the baseline model, which is linear and estimated\n"
681 "separately in each voxel (much as in 3dDeconvolve).\n"
682 "\n"
683 " -polort pnum = like in 3dDeconvolve [default is 'pnum' == 'A']\n"
684 "\n"
685 " -baseline bb = read 'bb' (a 1D file) for the baseline model;\n"
686 " this file can have 1 or more columns:\n"
687 " * if 1 column, then it is used in all slices\n"
688 " * if more than 1 column, then 'bb' must have\n"
689 " the same number of columns that the '-input'\n"
690 " dataset has slices, and each column in 'bb'\n"
691 " becomes part of the baseline model for only\n"
692 " the corresponding slice\n"
693 "\n"
694 "Response Model Arguments\n"
695 "------------------------\n"
696 "These arguments specify the response model. The Hemodynamic Response\n"
697 "Function (HRF) model is nonlinear. Given the HRF, the response to each\n"
698 "stimulus is modeled additively. Important notes:\n"
699 " * Each slice might have a different time offset, in which case the HRF\n"
700 " will be applied slightly differently in each slice.\n"
701 " * At least one '-stimtime' option must be given (or what would the\n"
702 " program be doing?).\n"
703 " * The same HRF applies to all voxels -- this is one distinction\n"
704 " between 3dNeocon and 3dDeconvolve, where the HRF is different\n"
705 " in each voxel. Only the amplitudes of the HRF fit vary between\n"
706 " voxels.\n"
707 " * The HRF itself has amplitude 1 (maximum absolute value). It is\n"
708 " the fit coefficients in each voxel that make the response model\n"
709 " vary in magnitude.\n"
710 " * Each time in the '-stimtime' inputs gets a separate amplitude estimate.\n"
711 " These need to be combined and/or contrasted in some other program,\n"
712 " such as 3dttest, to get a statistical map.\n"
713 "\n"
714 " -stimtime tname label HRFmodel\n"
715 "\n"
716 " 'tname' is the same format as 3dDeconvolve's '-stim_times' option\n"
717 " 'label' is a character string to identify the output coefficients\n"
718 " 'HRFmodel' specifies the form of the nonlinear response model; the\n"
719 " possibilities are\n"
720 " * 'GAMVAR' == the HRF has two parameters: 'tpeak' and 'fwhm';\n"
721 " HRF(t) = (t/bc)^b exp(b-t/c) * Heaviside(t)\n"
722 " b = (2.3*tpeak/fwhm)^2\n"
723 " c = tpeak/b\n"
724 " tpeak is allowed to range from 5 to 8 s;\n"
725 " fwhm is allowed to range from 4 to 8 s.\n"
726 " This HRF choice is appropriate for brief (under 2 s)\n"
727 " activations in response to each stimulus.\n"
728 " * 'IGAMVAR(d)' == similar to 'GAMVAR' but integrated over a duration\n"
729 " of 'd' seconds. This HRF choice is designed for block-design\n"
730 " FMRI experiments.\n"
731 " * 'CSPLIN(b,t,n)' == same as in 3dDeconvolve, with the caveat that\n"
732 " the maximum amplitude of the HRF will be one.\n"
733 " * 'DITTO' or 'ditto' == use the same model AND the same parameter set\n"
734 " as the previous '-stimtime' argument. In this way, the nonlinear\n"
735 " parameters for the HRFs for the two sets of time will be collapsed\n"
736 " into one set (e.g., both use the same value of 'tpeak' and 'fwhm').\n"
737 " Of course, the linear amplitudes for the two sets of times will be\n"
738 " separated.\n"
739 "\n"
740 " -threshtype hhh \n"
741 "\n"
742 " 'hhh' specifies the thresholding model used in the HRF analysis.\n"
743 " Given a set of HRF parameters, linear regression is used to fit\n"
744 " the baseline model and the baseline+response model in all voxels.\n"
745 " Those voxels that pass the threshold model (i.e., their response\n"
746 " fit is 'significant' in some sense) are used to judge the quality\n"
747 " of the overall fit. The choices for 'hhh' are\n"
748 " * 'RCLU(p,c)' == the nominal R^2-statistic (assuming white noise)\n"
749 " is computed in each voxel, and those voxels with a p-value at\n"
750 " or below 'p' are kept, if they are in a cluster of at least 'c'\n"
751 " such contiguous voxels. Subsequent to this thresholding, voxels\n"
752 " that fit better (have a larger R^2) are weighted more highly\n"
753 " in the HRF fitting objective function.\n"
754 " * At this time, there is no other thresholding model available.\n"
755 " * The default is currently 'RCLU(0.01,5)' which was just picked out\n"
756 " of thin air with no justification. This default is subject to\n"
757 " change!\n"
758 "\n"
759 "Output Arguments\n"
760 "----------------\n"
761 " -fitts fprefix = Output a fitted time series model dataset.\n"
762 " -errts eprefix = Output the residuals into a dataset.\n"
763 " -cbucket cprefix = Output the fit coefficients into a dataset.\n"
764 "\n"
765 " -verb = Increase the verbosity of the messages as the program runs.\n"
766 " -quiet = Turn off informational progress messages.\n"
767 ) ;
768
769 printf("\n*** AUTHOR = Zhark the Experimental, October 2007 ***\n\n" ) ;
770 return ;
771 }
772