1 #include "mrilib.h"
2
3 /*******************************************************************/
4 /********** 04 Mar 2003: Read a 1D file as an AFNI dataset *********/
5 /*******************************************************************/
6
7 /*-----------------------------------------------------------------*/
8 /*! Open a 1D file as an AFNI dataset. Each column is a separate
9 sub-brick. Ignore comments at this time. Data is not loaded
10 into bricks.
11 -------------------------------------------------------------------*/
12
THD_open_1D(char * pathname)13 THD_3dim_dataset * THD_open_1D( char *pathname )
14 {
15 THD_3dim_dataset *dset=NULL ;
16 MRI_IMAGE *flim ;
17 char prefix[1024] , *ppp , tname[12] , *pn ;
18 THD_ivec3 nxyz ;
19 THD_fvec3 dxyz , orgxyz ;
20 int nvals , lpn , flip=0 ;
21 FILE *fp ;
22
23 ENTRY("THD_open_1D") ;
24
25 /*-- open input file --*/
26
27 if( pathname == NULL || pathname[0] == '\0' ) RETURN(NULL);
28
29 /*-- check if it is a NIML-ish AFNI dataset;
30 if so, read it in that way instead of the 1D way [21 Mar 2003] --*/
31
32 if( strchr(pathname,'[') == NULL && strncmp(pathname,"1D:",3) != 0 &&
33 strchr(pathname,'{') == NULL && strchr(pathname,'\'') == NULL ){
34
35 pn = strdup(pathname) ; fp = fopen(pn,"r") ;
36
37 if( fp == NULL ){
38 char *p1 = strstr(pn,"1D") ; /* if can't open .1D, try .3D */
39 if( p1 != NULL ){
40 *p1 = '3' ; fp = fopen(pn,"r") ;
41 }
42 if( fp == NULL ){
43 fprintf(stderr,"** THD_open_1D(%s): can't open file\n",pathname);
44 free(pn); RETURN(NULL);
45 }
46 }
47 memset(prefix,0,133) ; fread(prefix,1,128,fp) ; fclose(fp) ;
48 if( strstr(prefix,"<AFNI_") != NULL && strstr(prefix,"dataset") != NULL ){
49 dset = THD_open_3D(pn) ;
50 if( dset != NULL && strcmp(pathname,pn) != 0 )
51 fprintf(stderr,"** THD_open_1D(%s): substituted %s\n",pathname,pn) ;
52 free(pn) ; return dset ;
53 }
54 }
55
56 /*-- otherwise, read it into an MRI_IMAGE, then mangle image into dataset --*/
57
58 lpn = strlen(pathname) ; pn = strdup(pathname) ;
59
60 #if 0
61 flip = (pn[lpn-1] == '\'') ; /* 12 Jul 2005: allow for tranposing input */
62 if( flip ) pn[lpn-1] = '\0' ;
63 #endif
64
65 flim = mri_read_1D(pn) ;
66 if( flim == NULL ){
67 fprintf(stderr,"** Can't read 1D dataset file %s\n",pn); free(pn); RETURN(NULL);
68 }
69 if( flip ){
70 MRI_IMAGE *qim = mri_transpose(flim); mri_free(flim); flim = qim;
71 }
72
73 /*-- make a dataset --*/
74
75 dset = EDIT_empty_copy(NULL) ; /* default (useless) dataset */
76
77 ppp = THD_trailname(pn,0) ; /* strip directory */
78 MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
79 ppp = strchr( prefix , '[' ) ;
80 if( ppp != NULL ) *ppp = '\0' ;
81
82 nxyz.ijk[0] = flim->nx ; dxyz.xyz[0] = 1.0 ; /* setup axes */
83 nxyz.ijk[1] = 1 ; dxyz.xyz[1] = 1.0 ;
84 nxyz.ijk[2] = 1 ; dxyz.xyz[2] = 1.0 ;
85
86 orgxyz.xyz[0] = 0.0 ; /* arbitrary origin */
87 orgxyz.xyz[1] = 0.0 ;
88 orgxyz.xyz[2] = 0.0 ;
89
90 nvals = flim->ny ; /* number of sub-bricks */
91
92 MCW_hash_idcode( pathname , dset ) ; /* 06 May 2005 */
93 dset->idcode.str[0] = 'A' ; /* overwrite 1st 3 bytes of IDcode */
94 dset->idcode.str[1] = '1' ;
95 dset->idcode.str[2] = 'D' ;
96
97 /* note we accept default orientation (RAI) here
98 (no use of ADN_xyzorient), since we actually
99 don't have any spatial structure to speak of */
100
101 EDIT_dset_items( dset ,
102 ADN_prefix , prefix ,
103 ADN_datum_all , MRI_float ,
104 ADN_nxyz , nxyz ,
105 ADN_xyzdel , dxyz ,
106 ADN_xyzorg , orgxyz ,
107 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
108 ADN_nvals , nvals ,
109 ADN_type , HEAD_ANAT_TYPE ,
110 ADN_func_type , ANAT_BUCK_TYPE ,
111 ADN_none ) ;
112
113 if( nvals > 1 && AFNI_yesenv("AFNI_1D_TIME") ){ /* pretend it is 3D+time */
114 char *eee = getenv("AFNI_1D_TIME_TR") ;
115 float dt = 0.0 ;
116 if( eee != NULL ) dt = strtod(eee,NULL) ;
117 if( dt <= 0.0 ) dt = 1.0 ;
118 EDIT_dset_items( dset ,
119 ADN_func_type, ANAT_EPI_TYPE ,
120 ADN_ntt , nvals ,
121 ADN_ttorg , 0.0 ,
122 ADN_ttdel , dt , /* TR */
123 ADN_ttdur , 0.0 ,
124 ADN_tunits , UNITS_SEC_TYPE ,
125 ADN_none ) ;
126 }
127
128 dset->dblk->diskptr->storage_mode = STORAGE_BY_1D ;
129 strcpy( dset->dblk->diskptr->brick_name , pathname ) ;
130
131 /*-- purge image data and return the empty dataset */
132
133 mri_free(flim) ; free(pn) ; RETURN(dset) ;
134 }
135
136 /*-----------------------------------------------------------------*/
137 /*! Load a 1D dataset's data into memory.
138 Called from THD_load_datablock() in thd_loaddblk.c.
139 -------------------------------------------------------------------*/
140
THD_load_1D(THD_datablock * dblk)141 void THD_load_1D( THD_datablock *dblk )
142 {
143 THD_diskptr *dkptr ;
144 MRI_IMAGE *flim ;
145 int nxyz , nbad,iv,nv ;
146 float *bar , *flar ;
147 char *pn ; int lpn , flip=0 ;
148
149 ENTRY("THD_load_1D") ;
150
151 /*-- open input [these errors should never occur] --*/
152
153 if( !ISVALID_DATABLOCK(dblk) ||
154 dblk->diskptr->storage_mode != STORAGE_BY_1D ||
155 dblk->brick == NULL ) EXRETURN ;
156
157 dkptr = dblk->diskptr ;
158 nxyz = dkptr->dimsizes[0] ;
159 nv = dkptr->nvals ;
160
161 if( nxyz*nv > 1000000 ) fprintf(stderr,"++ Reading %s\n",dkptr->brick_name) ;
162
163 pn = strdup(dkptr->brick_name) ; lpn = strlen(pn) ;
164 flip = (pn[lpn-1] == '\'') ; /* 12 Jul 2005: allow for tranposing input */
165 if( flip ) pn[lpn-1] = '\0' ;
166
167 flim = mri_read_1D(pn) ; free(pn) ;
168
169 if( flim == NULL ){
170 fprintf(stderr,"** THD_load_1D(%s): can't read file!\n",dkptr->brick_name) ;
171 EXRETURN ;
172 }
173 if( flip ){
174 MRI_IMAGE *qim = mri_transpose(flim); mri_free(flim); flim = qim;
175 }
176
177 /*-- allocate space for data --*/
178
179 if( nxyz != flim->nx || nv != flim->ny ){
180 fprintf(stderr,"** THD_load_1D(%s): nx or ny mismatch!\n",dkptr->brick_name) ;
181 fprintf(stderr,"** expect nx=%d; got nx=%d\n",nxyz,flim->nx) ;
182 fprintf(stderr,"** expect ny=%d; got ny=%d\n",nv ,flim->ny) ;
183 mri_free(flim) ; EXRETURN ;
184 }
185
186 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
187
188 /** malloc space for each brick separately **/
189
190 for( nbad=iv=0 ; iv < nv ; iv++ ){
191 if( DBLK_ARRAY(dblk,iv) == NULL ){
192 bar = AFMALL( float, DBLK_BRICK_BYTES(dblk,iv) ) ;
193 mri_fix_data_pointer( bar , DBLK_BRICK(dblk,iv) ) ;
194 if( bar == NULL ) nbad++ ;
195 }
196 }
197
198 /** if couldn't get them all, take our ball and go home in a snit **/
199
200 if( nbad > 0 ){
201 fprintf(stderr,"\n** failed to malloc %d 1D bricks out of %d\n\a",nbad,nv) ;
202 for( iv=0 ; iv < nv ; iv++ ){
203 if( DBLK_ARRAY(dblk,iv) != NULL ){
204 free(DBLK_ARRAY(dblk,iv)) ;
205 mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,iv) ) ;
206 }
207 }
208 mri_free(flim) ; EXRETURN ;
209 }
210
211 /** copy data from image to bricks **/
212
213 flar = MRI_FLOAT_PTR(flim) ;
214
215 for( iv=0 ; iv < nv ; iv++ ){
216 bar = DBLK_ARRAY(dblk,iv) ;
217 memcpy( bar , flar + iv*nxyz , sizeof(float)*nxyz ) ;
218 }
219
220 mri_free(flim) ; EXRETURN ;
221 }
222
223 /*------------------------------------------------------------------*/
224 /*! Write a dataset to disk as a 1D file.
225 Called from THD_write_3dim_dataset().
226 --------------------------------------------------------------------*/
227
THD_write_1D(char * sname,char * pname,THD_3dim_dataset * dset)228 void THD_write_1D( char *sname, char *pname , THD_3dim_dataset *dset )
229 {
230 char fname[THD_MAX_NAME] , *cpt ;
231 int iv,nv , nx,ny,nz,nxyz,ii,jj,kk , qq ;
232 FILE *fp=NULL ;
233 int binflag=0 ; char shp ; float val[1] ;
234 complex cval[1] ; int do_complex=0 ; char *dstr ;
235
236 ENTRY("THD_write_1D") ;
237
238 if( !ISVALID_DSET(dset) || !DSET_LOADED(dset) ) EXRETURN ;
239
240 nv = DSET_NVALS(dset) ; /* number of columns */
241 nx = DSET_NX(dset) ;
242 ny = DSET_NY(dset) ;
243 nz = DSET_NZ(dset) ; nxyz = nx*ny*nz ; /* number of rows */
244
245 /* make up a filename for output (into fname) */
246
247 cpt = DSET_PREFIX(dset) ;
248 if( (pname != NULL && *pname == '-') || /* 12 Jul 2005: to stdout */
249 (pname == NULL && (*cpt == '-' || strncmp(cpt,"stdout",6)==0)) ){
250 fp = stdout ; strcpy(fname,"stdout") ; binflag = 0 ;
251 }
252
253 /* 05 Mar 2008: special case to write as a 'pure' .1D file,
254 if filename ends in '.1D' or specifies stdout */
255
256 if( pname == NULL && AFNI_yesenv("AFNI_1D_TRANOUT") &&
257 (STRING_HAS_SUFFIX(cpt,".1D") || *cpt=='-' || strncmp(cpt,"stdout",6)==0) ){
258
259 MRI_IMAGE *qim = THD_dset_to_1Dmri(dset) ;
260 mri_write_1D(cpt,qim); mri_free(qim); EXRETURN ;
261 }
262
263 /* back to normal 3D mode of writing */
264
265 if( sname == NULL ) sname = DSET_DIRNAME(dset) ;
266
267 do_complex = fp != stdout &&
268 DSET_datum_constant(dset) && /* 24 Sep 2009: for LRF: */
269 DSET_BRICK_TYPE(dset,0) == MRI_complex ; /* allow complex output */
270
271 dstr = (do_complex) ? "complex" : "float" ;
272
273 if( fp == NULL ){
274 if( pname != NULL ){ /* have input prefix */
275 if( sname != NULL && *sname != '\0' ){ /* and input session (directory) */
276 strcpy(fname,sname) ;
277 ii = strlen(fname) ; if( fname[ii-1] != '/' ) strcat(fname,"/") ;
278 } else {
279 strcpy(fname,"./") ; /* don't have input session */
280 }
281 strcat(fname,pname) ;
282 } else { /* don't have input prefix */
283 if( sname != NULL && *sname != '\0' ){ /* use directory name -- 25 May 2009 */
284 strcpy(fname,sname) ;
285 ii = strlen(fname) ; if( fname[ii-1] != '/' ) strcat(fname,"/") ;
286 } else {
287 strcpy(fname,"./") ; /* don't have input session */
288 }
289 cpt = DSET_PREFIX(dset) ;
290 if( STRING_HAS_SUFFIX(cpt,".3D") || STRING_HAS_SUFFIX(cpt,".1D") )
291 strcat(fname,cpt) ;
292 else
293 strcpy(fname,DSET_BRIKNAME(dset)) ;
294
295 cpt = strchr(fname,'[') ;
296 if( cpt != NULL ) *cpt = '\0' ; /* without subscripts! */
297 }
298 ii = strlen(fname) ;
299 if( ii > 10 && strstr(fname,".BRIK") != NULL ){ /* delete +view.BRIK */
300 fname[ii-10] = '\0' ;
301 if( DSET_IS_1D(dset) || (ny==1 && nz==1) )
302 strcat(fname,".1D");
303 else
304 strcat(fname,".3D"); /* 21 Mar 2003 */
305 }
306
307 fp = fopen( fname , "w" ) ; if( fp == NULL ) EXRETURN ;
308 binflag = STRING_HAS_SUFFIX(fname,".3D") && AFNI_yesenv("AFNI_3D_BINARY") ;
309 }
310
311 strcpy( dset->dblk->diskptr->brick_name , fname ); /* 12 Jul 2005 */
312
313 /* are we going to write in binary? [03 Jun 2005] */
314
315 shp = (binflag) ? ' ' : '#' ;
316
317 /* write some dataset info as NIML-style header/comments */
318
319 if( fp != stdout )
320 fprintf(fp,
321 "%c <AFNI_3D_dataset\n"
322 "%c self_idcode = \"%s\"\n"
323 "%c ni_type = \"%d*%s\"\n" /* all columns are same type! */
324 "%c ni_dimen = \"%d,%d,%d\"\n"
325 "%c ni_delta = \"%g,%g,%g\"\n"
326 "%c ni_origin = \"%g,%g,%g\"\n"
327 "%c ni_axes = \"%s,%s,%s\"\n"
328 ,
329 shp ,
330 shp , dset->idcode.str ,
331 shp , nv , dstr ,
332 shp , nx,ny,nz ,
333 shp , DSET_DX(dset) , DSET_DY(dset) , DSET_DZ(dset) ,
334 shp , DSET_XORG(dset), DSET_YORG(dset), DSET_ZORG(dset),
335 shp , ORIENT_shortstr[dset->daxes->xxorient] ,
336 ORIENT_shortstr[dset->daxes->yyorient] ,
337 ORIENT_shortstr[dset->daxes->zzorient]
338 ) ;
339
340 if( fp != stdout && HAS_TIMEAXIS(dset) ){
341 float dt = DSET_TR(dset) ;
342 if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) dt *= 0.001 ;
343 fprintf(fp , "%c ni_timestep = \"%g\"\n" , shp,dt ) ;
344 }
345
346 if( fp != stdout && binflag )
347 fprintf(fp , " ni_form = \"binary.%s\"\n" ,
348 (NI_byteorder()==NI_LSB_FIRST) ? "lsbfirst" : "msbfirst" ) ;
349
350 /* do stataux for bricks, if any are present */
351
352 for( ii=iv=0 ; iv < nv ; iv++ )
353 if( DSET_BRICK_STATCODE(dset,iv) > 0 ) ii++ ;
354
355 if( fp != stdout && ii > 0 ){
356 fprintf(fp, "%c ni_stat = \"",shp) ;
357 for( iv=0 ; iv < nv ; iv++ ){
358 ii = DSET_BRICK_STATCODE(dset,iv) ;
359 if( ii <= 0 ){
360 fprintf(fp,"none") ;
361 } else {
362 fprintf(fp,"%s(",NI_stat_distname(ii)) ;
363 kk = NI_stat_numparam(ii) ;
364 for( jj=0 ; jj < kk ; jj++ ){
365 fprintf(fp,"%g",DSET_BRICK_STATPAR(dset,iv,jj)) ;
366 if( jj < kk-1 ) fprintf(fp,",") ;
367 }
368 fprintf(fp,")") ;
369 }
370 if( iv < nv-1 ) fprintf(fp,";") ;
371 }
372 fprintf(fp,"\"\n") ;
373 }
374
375 /* close NIML-style header */
376
377 if( fp != stdout ){
378 if( binflag ) fprintf(fp," >") ;
379 else fprintf(fp,"# >\n") ;
380 fflush(fp) ;
381 }
382
383 /* now write data */
384
385 for( ii=0 ; ii < nxyz ; ii++ ){ /* loop over voxels */
386
387 for( iv=0 ; iv < nv ; iv++ ){ /* loop over sub-bricks = columns */
388 switch( DSET_BRICK_TYPE(dset,iv) ){
389 default:
390 val[0] = 0.0f ;
391 break ;
392
393 case MRI_float:{
394 float *bar = DSET_ARRAY(dset,iv) ; val[0] = bar[ii] ;
395 }
396 break ;
397
398 case MRI_short:{
399 short *bar = DSET_ARRAY(dset,iv) ;
400 float fac = DSET_BRICK_FACTOR(dset,iv) ; if( fac == 0.0 ) fac = 1.0 ;
401 val[0] = fac*bar[ii] ;
402 }
403 break ;
404
405 case MRI_byte:{
406 byte *bar = DSET_ARRAY(dset,iv) ;
407 float fac = DSET_BRICK_FACTOR(dset,iv) ; if( fac == 0.0 ) fac = 1.0 ;
408 val[0] = fac*bar[ii] ;
409 }
410 break ;
411
412 /* below here, we convert complicated types to floats, losing data! */
413
414 case MRI_complex:{
415 complex *bar = DSET_ARRAY(dset,iv) ;
416 if( do_complex ) cval[0] = bar[ii] ;
417 else val[0] = CABS(bar[ii]) ;
418 }
419 break ;
420
421 case MRI_rgb:{
422 rgbyte *bar = DSET_ARRAY(dset,iv) ;
423 val[0] = (0.299*bar[ii].r+0.587*bar[ii].g+0.114*bar[ii].g) ;
424 }
425 break ;
426 } /* end of switch on sub-brick data type */
427
428 if( do_complex ){
429 if( binflag ) qq = fwrite( cval , sizeof(complex) , 1 , fp ) ;
430 else qq = fprintf( fp , " %g %g" , cval[0].r , cval[0].i ) ;
431 } else {
432 if( binflag ) qq = fwrite( val , sizeof(float) , 1 , fp ) ;
433 else qq = fprintf( fp , " %g" , val[0] ) ;
434 }
435
436 if( qq <= 0 ){ /* check for output error */
437 ERROR_message("THD_write_1D('%s') failure!",fname) ; goto DONE ;
438 }
439
440 } /* end of loop over sub-bricks */
441
442 if( !binflag ) fprintf(fp,"\n") ;
443
444 } /* end of loop over voxels */
445
446 /* NIML-style trailer */
447
448 DONE:
449 fflush(fp) ;
450
451 if( fp != stdout ){
452 fprintf(fp,"%c </AFNI_3D_dataset>\n",shp) ;
453 fclose(fp) ;
454 }
455
456 EXRETURN ;
457 }
458