1 #include "mrilib.h"
2
3 /*--------------------------------------------------------------------------*/
4
usage_3dExtractGroupInCorr(int detail)5 void usage_3dExtractGroupInCorr(int detail)
6 {
7 printf(
8 "Usage: 3dExtractGroupInCorr [options] AAA.grpincorr.niml\n"
9 "\n"
10 "This program breaks the collection of images from a GroupInCorr\n"
11 "file back into individual AFNI 3D+time datasets.\n"
12 "\n"
13 "Of course, only the data inside the mask used in 3dSetupGroupInCorr\n"
14 "is stored in the .data file, so only those portions of the input\n"
15 "files can be reconstructed :)\n"
16 "\n"
17 "The output datasets will be stored in float format, no matter what\n"
18 "the storage type of the original datasets or of the .data file.\n"
19 "\n"
20 "OPTION:\n"
21 "-------\n"
22 " -prefix PPP The actual dataset prefix with be the internal dataset\n"
23 " label with the string 'PPP_' pre-prended.\n"
24 " ++ Use NULL to skip the use of the prefix.\n"
25 "\n"
26 "Author -- RWCox -- May 2012\n"
27 ) ;
28 PRINT_COMPILE_DATE ;
29 return;
30 }
31
32 /*--------------------------------------------------------------------------*/
33
34 typedef struct {
35
36 int nvec ; /* number of vectors in a dataset */
37 int ndset ; /* number of datasets */
38 int *nvals ; /* nvals[i] = number of values in a vector in i-th dataset */
39 int nvals_max ; /* largest nvals[i] value */
40 int nvals_tot ; /* sum of nvals[i] values */
41 int datum ; /* 1 for sbyte, 2 for short */
42
43 int nuse , *use ; /* 07 Apr 2010: subset of datasets to use */
44
45 char *geometry_string ;
46 THD_3dim_dataset *tdset ; /* template dataset */
47 int nx,ny,nz,nvox ; /* number of nodes in template */
48
49 char *dfname ; /* data file name */
50 int *ivec ; /* ivec[i] = spatial index of i-th vector, i=0..nvec-1 */
51 float *fac ; /* fac[i] = scale factor for i-th dataset, i=0..ndset-1 */
52 short **sv ; /* sv[i] = short array [nvals[i]*nvec] for i-th dataset */
53 sbyte **bv ; /* bv[i] = sbyte array [nvals[i]*nvec] for i-th dataset */
54
55 char **dslab ; /* dslab[i] = label string for i-th dataset [23 May 2010] */
56
57 long long nbytes ; /* number of bytes in the data array */
58
59 /* Surface stuff ZSS Jan 09 2010 */
60
61 int nnode[2] ;
62 int ninmask[2] ;
63
64 } MRI_shindss ; /* short/sbyte indexed datasets */
65
66 /*--------------------------------------------------------------------------*/
67
68 #undef GQUIT
69 #define GQUIT(sss) \
70 do{ if( tdset != NULL ) DSET_delete(tdset) ; \
71 if( dfname != NULL ) free(dfname) ; \
72 if( geometry_string != NULL ) free(geometry_string) ; \
73 NI_free_element(nel) ; \
74 if( sss != NULL ) ERROR_message("EIC: file %s: %s",fname,(sss)) ; \
75 RETURN(NULL) ; \
76 } while(0)
77
78 /*--------------------------------------------------------------------------*/
79
80 static NI_element *nelshd = NULL ;
81
82 static const long long twogig = 2ll * 1024ll * 1024ll * 1024ll ; /* 2 GB */
83
84 /*----- read a PREFIX.grpincorr.niml file into a struct -----*/
85
GRINCOR_read_input(char * fname)86 MRI_shindss * GRINCOR_read_input( char *fname )
87 {
88 NI_element *nel=NULL ;
89 char *dfname=NULL , *atr ;
90 NI_float_array *facar ; NI_int_array *nvar, *nnode=NULL, *ninmask=NULL;
91 MRI_shindss *shd ;
92 long long nbytes_needed , nbytes_dfname=0 ; int fdes ;
93 void *var ; int ids , nvmax , nvtot ;
94 int datum , datum_size ;
95
96 char *geometry_string=NULL ;
97 THD_3dim_dataset *tdset=NULL; int nvox;
98 int no_ivec=0 , *ivec=NULL , *nvals=NULL , nvec,ndset ; float *fac=NULL ;
99 NI_str_array *slabar=NULL ;
100
101 if( fname == NULL || *fname == '\0' ) GQUIT(NULL) ;
102
103 ENTRY("GRINCOR_read_input") ;
104
105 /* get data element */
106
107 if (!THD_is_ondisk(fname))
108 GQUIT("not on disk") ;
109
110 nelshd = nel = NI_read_element_fromfile(fname) ;
111
112 if( nel == NULL || nel->type != NI_ELEMENT_TYPE )
113 GQUIT("not properly formatted") ;
114 if( strcmp(nel->name,"3dGroupInCorr") != 0 )
115 GQUIT("data element name is not '3dGroupInCorr'") ;
116
117 /* no data vector ==> using all voxels */
118
119 no_ivec = ( nel->vec_num < 1 ||
120 nel->vec_len < 1 || nel->vec_typ[0] != NI_INT ) ;
121
122 /* number of vectors in each dataset */
123
124 atr = NI_get_attribute(nel,"nvec");
125 if( atr == NULL ) GQUIT("nvec attribute missing?") ;
126 nvec = (int)strtod(atr,NULL) ;
127 if( nvec < 2 || (!no_ivec && nel->vec_len != nvec) )
128 GQUIT("nvec attribute has illegal value") ;
129
130 /* number of datasets */
131
132 atr = NI_get_attribute(nel,"ndset");
133 if( atr == NULL ) GQUIT("ndset attribute missing") ;
134 ndset = (int)strtod(atr,NULL) ;
135 if( ndset < 1 ) GQUIT("ndset attribute has illegal value") ;
136
137 /* number of time points in each dataset (varies with dataset) */
138
139 atr = NI_get_attribute(nel,"nvals");
140 if( atr == NULL ) GQUIT("nvals attribute missing") ;
141 nvar = NI_decode_int_list(atr,",") ;
142 if( nvar == NULL || nvar->num < ndset )
143 GQUIT("nvals attribute doesn't match ndset") ;
144 nvals = nvar->ar ; nvar->ar = NULL ; NI_delete_int_array(nvar) ;
145
146 nvmax = nvtot = nvals[0] ;
147 for( ids=1 ; ids < ndset ; ids++ ){ /* Feb 2011 */
148 nvtot += nvals[ids] ;
149 if( nvals[ids] > nvmax ) nvmax = nvals[ids] ;
150 }
151
152 /* dataset labels [23 May 2010] */
153
154 atr = NI_get_attribute(nel,"dset_labels") ;
155 if( atr != NULL ){
156 slabar = NI_decode_string_list(atr,";,") ;
157 if( slabar == NULL || slabar->num < ndset )
158 GQUIT("dset_labels attribute invalid") ;
159 }
160
161 /* datum of datasets */
162
163 atr = NI_get_attribute(nel,"datum") ;
164 if( atr != NULL && strcasecmp(atr,"byte") == 0 ){
165 datum = 1 ; datum_size = sizeof(sbyte) ;
166 } else {
167 datum = 2 ; datum_size = sizeof(short) ;
168 }
169
170 /* number of bytes needed:
171 sizeof(datum) * number of vectors per dataset
172 * number of datasets
173 * sum of per dataset vector lengths */
174
175 nbytes_needed = 0 ;
176 for( ids=0 ; ids < ndset ; ids++ ) nbytes_needed += nvals[ids] ;
177 nbytes_needed *= ((long long)nvec) * datum_size ;
178
179 if( nbytes_needed >= twogig &&
180 ( sizeof(void *) < 8 || sizeof(size_t) < 8 ) ) /* too much for 32-bit */
181 GQUIT("datafile size exceeds 2 GB -- you need a 64-bit computer!") ;
182
183 /* scale factor for each dataset */
184
185 atr = NI_get_attribute(nel,"fac") ;
186 if( atr == NULL ) GQUIT("fac attribute missing") ;
187 facar = NI_decode_float_list(atr,",") ;
188 if( facar == NULL || facar->num < ndset )
189 GQUIT("can't decode fac attribute") ;
190 fac = facar->ar ; facar->ar = NULL ; NI_delete_float_array(facar) ;
191
192 for( ids=0 ; ids < ndset ; ids++ ) if( fac[ids] <= 0.0f ) fac[ids] = 1.0f ;
193
194 /* grid definition */
195
196 atr = NI_get_attribute(nel,"geometry") ;
197 if( atr == NULL ) GQUIT("geometry attribute missing") ;
198 geometry_string = strdup(atr) ;
199 tdset = EDIT_geometry_constructor( geometry_string , "GrpInCorr" ) ;
200 if( tdset == NULL ) GQUIT("can't decode geometry attribute") ;
201 nvox = DSET_NVOX(tdset) ;
202 if( no_ivec && nvox != nvec )
203 GQUIT("geometry attribute doesn't match nvec attribute") ;
204 if( !no_ivec && nvox < nvec )
205 GQUIT("geometry attribute specifies too few voxels") ;
206
207 /* name of data file: check its size against what's needed */
208
209 #if 0
210 atr = NI_get_attribute(nel,"datafile") ;
211 if( atr != NULL ){
212 dfname = strdup(atr) ; nbytes_dfname = THD_filesize(dfname) ;
213 if( nbytes_dfname <= 0 && strstr(dfname,"/") != NULL ){
214 char *tnam = THD_trailname(atr,0) ;
215 nbytes_dfname = THD_filesize(tnam) ;
216 if( nbytes_dfname > 0 ){ free(dfname); dfname = strdup(tnam); }
217 }
218 }
219 #endif
220 if( nbytes_dfname <= 0 && strstr(fname,".niml") != NULL ){
221 if( dfname != NULL ) free(dfname) ;
222 dfname = strdup(fname) ; strcpy(dfname+strlen(dfname)-5,".data") ;
223 nbytes_dfname = THD_filesize(dfname) ;
224 }
225 if( nbytes_dfname <= 0 ){
226 char mess[THD_MAX_NAME+256] ;
227 sprintf(mess,"datafile is missing (%s)",dfname) ; GQUIT(mess) ;
228 } else if( nbytes_dfname < nbytes_needed ){
229 char mess[THD_MAX_NAME+1024] ;
230 sprintf(mess,"datafile %s has %s bytes but needs at least %s",
231 dfname ,
232 commaized_integer_string(nbytes_dfname) ,
233 commaized_integer_string(nbytes_needed) ) ;
234 GQUIT(mess) ;
235 } else {
236 INFO_message("EIC: data file %s found with %s bytes of data",
237 dfname , commaized_integer_string(nbytes_dfname) ) ;
238 }
239 fdes = open( dfname , O_RDWR ) ;
240 if( fdes < 0 ){
241 char mess[THD_MAX_NAME+256] ;
242 sprintf(mess,"can't open datafile (%s)",dfname) ; GQUIT(mess) ;
243 }
244 NI_set_attribute( nelshd , "datafile" , dfname ) ;
245
246 /* ivec[i] is the voxel spatial index of the i-th vector */
247
248 if( no_ivec ){
249 ivec = NULL ; /* means all voxels: ivec[i] == i */
250 } else {
251 ivec = (int *)nel->vec[0] ; /* copy pointer */
252 nel->vec[0] = NULL ; /* NULL out in element so won't be free-ed */
253 }
254
255 /* And stuff for LR surface pairs ZSS Jan 09*/
256 if ((atr=NI_get_attribute(nel,"LRpair_nnode"))) {
257 nnode = NI_decode_int_list(atr,",") ;
258 }
259 if ((atr=NI_get_attribute(nel,"LRpair_ninmask"))) {
260 ninmask = NI_decode_int_list(atr,",") ;
261 }
262
263 /* create output struct */
264
265 shd = (MRI_shindss *)malloc(sizeof(MRI_shindss)) ;
266
267 shd->nvals = nvals ; shd->nvals_max = nvmax ; shd->nvals_tot = nvtot ;
268 shd->nvec = nvec ;
269 shd->ndset = ndset ;
270
271 shd->geometry_string = geometry_string ;
272 shd->tdset = tdset ;
273 shd->dfname = dfname ;
274 shd->nvox = nvox ;
275 shd->nx = DSET_NX(tdset); shd->ny = DSET_NY(tdset); shd->nz = DSET_NZ(tdset);
276
277 shd->ivec = ivec ;
278 shd->fac = fac ;
279
280 /* and surface fields... ZSS Jan 09 */
281 if (nnode) {
282 if (nnode->num != 2) GQUIT("LRpair_nnode must have 2 values");
283 shd->nnode[0] = nnode->ar[0];
284 shd->nnode[1] = nnode->ar[1];
285 NI_delete_int_array(nnode); nnode=NULL;
286 } else {
287 shd->nnode[0] = shd->nnode[1] = -1 ;
288 }
289 if (ninmask) {
290 if (ninmask->num != 2) GQUIT("LRpair_ninmask must have 2 values");
291 shd->ninmask[0] = ninmask->ar[0];
292 shd->ninmask[1] = ninmask->ar[1];
293 NI_delete_int_array(ninmask); ninmask=NULL;
294 } else {
295 shd->ninmask[0] = shd->ninmask[1] = -1 ;
296 }
297
298 /*--- 07 Apr 2010: setup default use list (all of them) ---*/
299
300 shd->nuse = ndset ;
301 shd->use = (int *)malloc(sizeof(int)*ndset) ;
302 for( ids=0 ; ids < ndset ; ids++ ) shd->use[ids] = ids ;
303
304 shd->dslab = (slabar != NULL) ? slabar->str : NULL ; /* 23 May 2010 */
305
306 /*--- now have to map data from disk ---*/
307
308 var = mmap( 0 , (size_t)nbytes_needed ,
309 PROT_WRITE , THD_MMAP_FLAG , fdes , 0 ) ;
310 close(fdes) ; /* close file descriptor does not unmap data */
311
312 if( var == (void *)(-1) ){ /* this is bad */
313 ERROR_message(
314 "EIC: file %s: can't mmap() datafile -- memory space exhausted?" , dfname ) ;
315 free(shd) ; RETURN(NULL) ;
316 }
317
318 /*-- create array of pointers to each dataset's data array --*/
319
320 shd->datum = datum ;
321
322 if( datum == 2 ){ /* shorts */
323 shd->sv = (short **)malloc(sizeof(short *)*ndset) ;
324 shd->bv = NULL ;
325 shd->sv[0] = (short *)var ;
326 for( ids=1 ; ids < ndset ; ids++ )
327 shd->sv[ids] = shd->sv[ids-1] + nvals[ids-1]*nvec ;
328 } else { /* sbytes */
329 shd->sv = NULL ;
330 shd->bv = (sbyte **)malloc(sizeof(sbyte *)*ndset) ;
331 shd->bv[0] = (sbyte *)var ;
332 for( ids=1 ; ids < ndset ; ids++ )
333 shd->bv[ids] = shd->bv[ids-1] + nvals[ids-1]*nvec ;
334 }
335
336 shd->nbytes = nbytes_needed ;
337 RETURN(shd) ;
338 }
339
340 #undef GQUIT
341
342 /*----------------------------------------------------------------------------*/
343 /* Create a vectim struct from 1 dataset inside the collection */
344
GRINCOR_extract_vectim_short(MRI_shindss * shd,int ids)345 MRI_vectim * GRINCOR_extract_vectim_short( MRI_shindss *shd , int ids )
346 {
347 MRI_vectim *mv ;
348 long long nvec=shd->nvec , nvals=shd->nvals[ids] , ii,nvv ;
349 float fac=shd->fac[ids] , *fv ;
350 short *sv = shd->sv[ids] ;
351
352 ENTRY("GRINCOR_extract_vectim_short") ;
353
354 MAKE_VECTIM( mv , nvec , nvals ) ;
355 fv = mv->fvec ;
356 nvv = nvec * nvals ;
357 for( ii=0 ; ii < nvv ; ii++ ) fv[ii] = fac * sv[ii] ;
358 RETURN(mv) ;
359 }
360
GRINCOR_extract_vectim_sbyte(MRI_shindss * shd,int ids)361 MRI_vectim * GRINCOR_extract_vectim_sbyte( MRI_shindss *shd , int ids )
362 {
363 MRI_vectim *mv ;
364 long long nvec=shd->nvec , nvals=shd->nvals[ids] , ii,nvv ;
365 float fac=shd->fac[ids] , *fv ;
366 sbyte *sv = shd->bv[ids] ;
367
368 ENTRY("GRINCOR_extract_vectim_sbyte") ;
369
370 MAKE_VECTIM( mv , nvec , nvals ) ;
371 fv = mv->fvec ;
372 nvv = nvec * nvals ;
373 for( ii=0 ; ii < nvv ; ii++ ) fv[ii] = fac * sv[ii] ;
374 RETURN(mv) ;
375 }
376
GRINCOR_extract_vectim(MRI_shindss * shd,int ids)377 MRI_vectim * GRINCOR_extract_vectim( MRI_shindss *shd , int ids )
378 {
379 MRI_vectim *mv ;
380 if( shd->datum == 1 ) mv = GRINCOR_extract_vectim_sbyte(shd,ids) ;
381 else mv = GRINCOR_extract_vectim_short(shd,ids) ;
382 return mv ;
383 }
384
385 /*----------------------------------------------------------------------------*/
386
GRINCOR_extract_dataset(MRI_shindss * shd,int ids,char * pref)387 THD_3dim_dataset * GRINCOR_extract_dataset( MRI_shindss *shd, int ids, char *pref )
388 {
389 MRI_vectim *mv ;
390 THD_3dim_dataset *dset ;
391 char prefix[THD_MAX_NAME] ;
392 int iv , nvals=shd->nvals[ids] ;
393 static int nds=0 ;
394
395 ENTRY("GRINCOR_extract_dataset") ;
396
397 STATUS("extract vectim") ;
398 mv = GRINCOR_extract_vectim( shd , ids ) ;
399
400 STATUS("create empty copy of template") ;
401 dset = EDIT_empty_copy( shd->tdset ) ;
402
403 STATUS("edit prefix") ;
404 prefix[0] = '\0' ;
405 if( pref != NULL && *pref != '\0' ){ strcpy(prefix,pref) ; strcat(prefix,"_") ; }
406 if( shd->dslab != NULL && shd->dslab[ids] != NULL ){
407 strcat(prefix,shd->dslab[ids]) ;
408 } else {
409 nds++ ;
410 sprintf(prefix+strlen(prefix),"%03d",nds) ;
411 }
412
413 STATUS("edit empty copy header") ;
414 EDIT_dset_items( dset ,
415 ADN_prefix , prefix ,
416 ADN_nvals , nvals ,
417 ADN_ntt , nvals ,
418 ADN_ttdel , 1.0 ,
419 ADN_tunits , UNITS_SEC_TYPE ,
420 ADN_brick_fac , NULL ,
421 ADN_type , HEAD_FUNC_TYPE ,
422 ADN_func_type , FUNC_FIM_TYPE ,
423 ADN_none ) ;
424
425 STATUS("create empty float bricks") ;
426 for( iv=0 ; iv < nvals ; iv++ )
427 EDIT_substitute_brick( dset , iv , MRI_float , NULL ) ;
428
429 STATUS("copy index vector") ;
430 if( shd->ivec != NULL ){
431 memcpy( mv->ivec , shd->ivec , sizeof(int)*shd->nvec ) ;
432 } else {
433 for( iv=0 ; iv < shd->nvec ; iv++ ) mv->ivec[iv] = iv ;
434 }
435
436 STATUS("convert vectim to dset") ;
437 THD_vectim_to_dset( mv , dset ) ;
438
439 STATUS("destroy vectim") ;
440 VECTIM_destroy( mv ) ;
441
442 RETURN(dset) ;
443 }
444
445 /*----------------------------------------------------------------------------*/
446
main(int argc,char * argv[])447 int main( int argc , char *argv[] )
448 {
449 MRI_shindss *shd ;
450 int ids , nopt , kk ;
451 char *prefix = "EIC" ;
452 char *fname=NULL , *buf ;
453 THD_3dim_dataset *dset ;
454
455 /*--- official AFNI startup stuff ---*/
456
457 mainENTRY("3dExtractGroupInCorr"); machdep();
458 AFNI_logger("3dExtractGroupInCorr",argc,argv);
459 PRINT_VERSION("3dExtractGroupInCorr"); AUTHOR("RW Cox");
460
461 /*-- process options --*/
462
463 nopt = 1 ;
464 while( nopt < argc && argv[nopt][0] == '-' ){
465
466 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
467 nopt++ ;
468 if( strcasecmp(argv[nopt],"NULL") == 0 ) prefix = NULL ;
469 else prefix = strdup(argv[nopt]) ;
470 nopt++ ; continue ;
471 }
472
473 if( strcmp(argv[nopt],"-help") == 0 || strcmp(argv[nopt],"-h") == 0 ){
474 usage_3dExtractGroupInCorr(2) ;
475 exit(0) ;
476 }
477
478 ERROR_message("Unknown option: '%s'",argv[nopt]) ;
479 suggest_best_prog_option(argv[0], argv[nopt]);
480 exit(1);
481 }
482
483 if( argc < 2 ){ usage_3dExtractGroupInCorr(2) ; exit(0) ; }
484
485 /* check for errors */
486
487 if( nopt >= argc ) ERROR_exit("No input filename on command line?!") ;
488
489 /*-- read input file --*/
490
491 fname = strdup(argv[nopt]) ;
492 if( STRING_HAS_SUFFIX(fname,".data") ){
493 strcpy(fname+strlen(fname)-5,".niml") ;
494 WARNING_message("EIC: Replaced '.data' with '.niml' in filename") ;
495 } else if( STRING_HAS_SUFFIX(fname,".grpincorr") ){
496 fname = (char *)realloc(fname,strlen(fname)+16) ;
497 strcat(fname,".niml") ;
498 INFO_message("EIC: Added '.niml' to end of filename") ;
499 } else if( STRING_HAS_SUFFIX(fname,".grpincorr.") ){
500 fname = (char *)realloc(fname,strlen(fname)+16) ;
501 strcat(fname,"niml") ;
502 INFO_message("EIC: Added 'niml' to end of filename") ;
503 }
504 shd = GRINCOR_read_input( fname ) ;
505 if( shd == NULL ) ERROR_exit("EIC: Cannot continue after input error") ;
506 INFO_message("EIC: file opened, contains %d datasets, %d time series, %s bytes",
507 shd->ndset , shd->nvec , commaized_integer_string(shd->nbytes) ) ;
508
509 /*-- process input file --*/
510
511 fprintf(stderr,"++ %d datasets: ",shd->ndset) ;
512 for( ids=0 ; ids < shd->ndset ; ids++ ){
513 fprintf(stderr,"%d",ids+1) ;
514 dset = GRINCOR_extract_dataset( shd, ids, prefix ) ; fprintf(stderr,".") ;
515 DSET_write(dset) ;
516 DSET_delete(dset) ;
517 }
518 fprintf(stderr,"\n") ; exit(0) ;
519 }
520