1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*! \file
8 This file contains all the functions for reading image files.
9 It is primarily used by to3d.c. It reads 2D images (into MRI_IMAGE struct)
10 and arrays of 2D images (into MRI_IMARR struct).
11 */
12
13 /*--------------------------------------------------------*/
14 /*** 7D SAFE (but most routines only return 2D images!) ***/
15 /*--------------------------------------------------------*/
16
17 #include <stdio.h>
18 #include <ctype.h>
19 #include <string.h>
20 #include <sys/stat.h>
21
22 /*** for non ANSI compilers ***/
23
24 #ifndef SEEK_SET
25 #define SEEK_SET 0
26 #endif
27
28 #include "mrilib.h"
29 #include "ge4_header.h"
30
31 #undef NLL
32 #define NLL 32222 /* lbuf length below -- should be enuf */
33
34 /*---------------------------------------------------------------*/
35 static MRI_IMAGE * mri_try_mri( FILE * , int * ) ; /* prototypes */
36 static MRI_IMAGE * mri_try_7D ( FILE * , int * ) ;
37 static MRI_IMAGE * mri_try_pgm( FILE * , int * ) ;
38 static int check_dicom_magic_num( char * ) ;
39 /*---------------------------------------------------------------*/
40
41 /*! Global variable to signal image orientation, if possible. */
42
43 char MRILIB_orients[8] = "\0\0\0\0\0\0\0\0" ; /* 12 Mar 2001 */
44
45 /*! Global variable to signal image slice offset, if possible. */
46
47 float MRILIB_zoff = 0.0 ;
48
49 /*! Global variable to signal image TR, if possible. */
50
51 float MRILIB_tr = 0.0 ; /* 03 Dec 2001 */
52
53 /*! Global variable to signal image x offset, if possible. */
54
55 float MRILIB_xoff = 0.0 ; /* 07 Dec 2001 */
56
57 /*! Global variable to signal image y offset, if possible. */
58
59 float MRILIB_yoff = 0.0 ;
60
61 /*! Global variable saying whether to use MRILIB_xoff. */
62
63 int use_MRILIB_xoff = 0 ;
64
65 /*! Global variable saying whether to use MRILIB_yoff. */
66
67 int use_MRILIB_yoff = 0 ;
68
69 /*! Global variable saying whether to use MRILIB_zoff. */
70
71 int use_MRILIB_zoff = 0 ;
72
73 /*! Global variable saying whether to use MRILIB_xcos. */
74
75 int use_MRILIB_xcos = 0 ;
76
77 /*! Global vector pointing in direction of x-axis. */
78
79 float MRILIB_xcos[3] = { 1.0 , 0.0 , 0.0 } ;
80
81 /*! Global variable saying whether to use MRILIB_ycos. */
82
83 int use_MRILIB_ycos = 0 ;
84
85 /*! Global vector pointing in direction of y-axis. */
86
87 float MRILIB_ycos[3] = { 0.0 , 1.0 , 0.0 } ;
88
89 /*! Global variable saying whether to use MRILIB_zcos. */
90
91 int use_MRILIB_zcos = 0 ;
92
93 /*! Global vector pointing in direction of z-axis. */
94
95 float MRILIB_zcos[3] = { 0.0 , 0.0 , 1.0 } ;
96
97 /*! Global variable saying whether to use MRILIB_slicespacing. */
98
99 int use_MRILIB_slicespacing = 0 ;
100
101 /*! Global variable giving the spacing between slice centers. */
102
103 float MRILIB_slicespacing = 0.0 ;
104
105 /*! Global variable saying whether to use DICOM matrix below. */
106
107 int use_MRILIB_dicom_matrix = 0 ;
108
109 /*! Global variable defining 3D image position and orientation. */
110
111 mat44 MRILIB_dicom_matrix ;
112
113 /*! Global variable counting number of DICOM files read. */
114
115 int MRILIB_dicom_count = 0 ;
116
117 /*! Global variable reporting signed 16-bit overflow (into unsigned). */
118
119 int MRILIB_dicom_s16_overflow = 0 ;
120
121 /*! Global variable containing DICOM header info for each file. */
122
123 AFD_dicom_header **MRILIB_dicom_header = NULL ;
124
125 /*-----------------------------------------------------------------*/
126
127 /*! Database of preset file sizes for auto-use of 3D.
128
129 If( head < 0 ) then file length must match size
130 else file length must == n*size + head for some
131 integer n, which will be the number of slices
132 to read from the file.
133 */
134
135 typedef struct {
136 int size ; /*!< file size in bytes if head < 0; image size in bytes if head >= 0 */
137 int head ; /*!< global header size if >= 0 */
138 char * prefix ; /*!< character string to prefix to filename */
139 } MCW_imsize ;
140
141 /*! Max number of preset file sizes to allow. */
142
143 #define MAX_MCW_IMSIZE 99
144
145 /*! Array of preset file sizes to use when reading image files. */
146
147 static MCW_imsize imsize[MAX_MCW_IMSIZE] ;
148
149 /*! If this < 0 ==> must initialize array of preset file sizes. */
150
151 static int MCW_imsize_good = -1 ;
152
153 /*------------------------------------------------------------------------*/
154 /* Siemens slice timing info and functions 13 Apr 2011 [rickr] */
155 /* note: nalloc is not propagated, nused == nalloc here */
156
157 int g_siemens_timing_nused = 0;
158 float * g_siemens_timing_times = NULL;
159 int g_siemens_timing_units = UNITS_MSEC_TYPE; /* not yet varying */
160
cleanup_g_siemens_times()161 static int cleanup_g_siemens_times()
162 {
163 if( ! g_siemens_timing_nused ) return 0;
164
165 if( g_siemens_timing_times ) free(g_siemens_timing_times);
166 g_siemens_timing_nused = 0;
167 g_siemens_timing_times = NULL;
168
169 return 0;
170 }
171
alloc_g_siemens_times(int ntimes)172 static int alloc_g_siemens_times(int ntimes)
173 {
174 if( ntimes <= 0 ) return cleanup_g_siemens_times();
175
176 /* if there is nothing to allocate, we're good */
177 if( ntimes == g_siemens_timing_nused ) return 0;
178
179 /* actually allocate memory */
180 g_siemens_timing_times = (float *)realloc(g_siemens_timing_times,
181 ntimes * sizeof(float));
182 if( ! g_siemens_timing_times ) {
183 fprintf(stderr,"** siemens AGST: failed to alloc %d floats\n", ntimes);
184 cleanup_g_siemens_times();
185 return 1;
186 }
187
188 /* and note the new size */
189 g_siemens_timing_nused = ntimes;
190
191 return 0;
192 }
193
populate_g_siemens_times(int tunits)194 int populate_g_siemens_times(int tunits)
195 {
196 float * d_times, tfac = 1.0;
197 int d_nalloc, d_nused; /* to be copied from mri_dicom_hdr.c */
198 int index;
199
200 ENTRY("populate_g_siemens_times");
201
202 if( mri_siemens_slice_times(&d_nalloc, &d_nused, &d_times) ) {
203 /* odd failure */
204 fprintf(stderr,"** PGST: odd failure getting siemens slice times\n");
205 cleanup_g_siemens_times();
206 RETURN(1);
207 }
208
209 /* allocate any needed memory */
210 if( alloc_g_siemens_times(d_nused) ) RETURN(1);
211
212 /* if no times, we're done */
213 if( d_nused == 0 ) RETURN(0);
214
215 /* siemens units are ms, so if we want seconds, divide by 1000 */
216 if( tunits == UNITS_SEC_TYPE ) tfac = 0.001;
217 else if ( tunits == UNITS_MSEC_TYPE ) tfac = 1.0;
218 else fprintf(stderr,"** PGST: bad time units %d\n", tunits);
219
220 /* copy the data (one by one, to allow for time scalar) */
221 for( index = 0; index < d_nused; index++ )
222 g_siemens_timing_times[index] = d_times[index] * tfac;
223
224 RETURN(0);
225 }
226
227 /* given the number of slices, the volume TR and the time units:
228 * 1. verify that nz matches the number of times
229 * 2. verify that 0 <= min, max <= TR
230 * 3. maybe output which order the times seem to be in (future?)
231 *
232 * if range_check, times outside TR are an error
233 *
234 * if verb, output any errors
235 * if verb > 1, output time pattern to screen
236 *
237 * return 1 if valid, 0 otherwise
238 */
valid_g_siemens_times(int nz,float TR,int range_check,int verb)239 int valid_g_siemens_times(int nz, float TR, int range_check, int verb)
240 {
241 float min, max, * times = g_siemens_timing_times;
242 int ind, decimals=3;
243
244 ENTRY("test_g_siemens_times");
245
246 if( nz != g_siemens_timing_nused ) {
247 if(verb) {
248 fprintf(stderr,"** ERROR: have %d siemens times but %d slices\n",
249 g_siemens_timing_nused, nz);
250 fprintf(stderr," Consider 'dicom_hdr -slice_times' for details.\n");
251 }
252 RETURN(0);
253 }
254
255 if( nz < 1 ) RETURN(1);
256
257 /* get min and max */
258 min = max = times[0];
259 for( ind = 1; ind < nz; ind++ ) {
260 if( times[ind] < min ) min = times[ind];
261 if( times[ind] > max ) max = times[ind];
262 }
263
264 if( verb > 1 ) { /* print the times */
265 if( max > 100 ) decimals = 1;
266 else decimals = 3;
267 printf("-- using Siemens slice timing (%d) :", nz);
268 for( ind = 0; ind < nz; ind++ ) printf(" %.*f", decimals, times[ind]);
269 putchar('\n');
270 }
271
272 /* use stdout to report issues here, to stick with times report */
273 if( min < 0.0 ) {
274 if(verb) printf("** min slice time %.*f outside TR range [0.0, %.*f]\n",
275 decimals, min, decimals, TR);
276 if( range_check ) RETURN(0); /* bad times, man 9 Nov 2016 [rickr] */
277 }
278 else if( max > TR ) {
279 if(verb) printf("** max slice time %.*f outside TR range [0.0, %.*f]\n",
280 decimals, max, decimals, TR);
281 if( range_check ) RETURN(0); /* bad times, man */
282 }
283
284 RETURN(1);
285 }
286
get_and_display_siemens_times(void)287 int get_and_display_siemens_times(void)
288 {
289 float * times;
290 int ind, ntimes;
291
292 ENTRY("get_and_display_siemens_times");
293
294 if( populate_g_siemens_times(UNITS_MSEC_TYPE) ) RETURN(1);
295
296 ntimes = g_siemens_timing_nused;
297 times = g_siemens_timing_times;
298
299 if( ntimes <= 0 ) {
300 printf("-- no Siemens timing found\n");
301 RETURN(0);
302 }
303
304 printf("-- Siemens timing (%d entries):", ntimes);
305 for( ind = 0; ind < ntimes; ind++ ) printf(" %.1f", times[ind]);
306 putchar('\n');
307
308 RETURN(0);
309 }
310
311 /* end siemens slice timing globals and functions */
312 /*------------------------------------------------------------------------*/
313
314
315 /*---------------------------------------------------------------*/
316
317 #undef swap_4
318 #undef swap_2
319
320 /*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
321
swap_4(void * ppp)322 static void swap_4(void *ppp)
323 {
324 unsigned char *pntr = (unsigned char *) ppp ;
325 unsigned char b0, b1, b2, b3;
326
327 b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
328 *pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
329 }
330
331 /*---------------------------------------------------------------*/
332
333 /*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
334
swap_8(void * ppp)335 static void swap_8(void *ppp)
336 {
337 unsigned char *pntr = (unsigned char *) ppp ;
338 unsigned char b0, b1, b2, b3;
339 unsigned char b4, b5, b6, b7;
340
341 b0 = *pntr ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
342 b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
343
344 *pntr = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
345 *(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
346 }
347
348 /*---------------------------------------------------------------*/
349
350 /*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
351
swap_2(void * ppp)352 static void swap_2(void *ppp)
353 {
354 unsigned char *pntr = (unsigned char *) ppp ;
355 unsigned char b0, b1;
356
357 b0 = *pntr; b1 = *(pntr+1);
358 *pntr = b1; *(pntr+1) = b0;
359 }
360
361 /******************************************************************/
362
363 /*! Earliest image reading function in the AFNI package.
364 Reads a single 2D image.
365
366 \param fname is the name of the file to try to read
367 \return NULL if an image couldn't be read, otherwise
368 a pointer to an MRI_IMAGE with data, dimensions, etc.
369 */
370
mri_read(char * fname)371 MRI_IMAGE *mri_read( char *fname )
372 {
373 FILE *imfile ;
374 MRI_IMAGE *im=NULL ;
375 int length , skip=0 , swap=0 ;
376 void *data ;
377
378 ENTRY("mri_read") ;
379
380 if( fname == NULL || *fname == '\0' ) RETURN(NULL) ; /* bad user */
381
382 /**-- 27 Apr 2005: check here for special filenames --**/
383
384 if( strstr(fname,".jpg" ) != NULL || /* various formats */
385 strstr(fname,".JPG" ) != NULL || /* that we convert */
386 strstr(fname,".jpeg") != NULL || /* to PPG/PGM using */
387 strstr(fname,".JPEG") != NULL || /* external filters */
388 strstr(fname,".gif" ) != NULL ||
389 strstr(fname,".GIF" ) != NULL ||
390 strstr(fname,".tif" ) != NULL ||
391 strstr(fname,".TIF" ) != NULL ||
392 strstr(fname,".tiff") != NULL ||
393 strstr(fname,".TIFF") != NULL ||
394 strstr(fname,".heic") != NULL ||
395 strstr(fname,".HEIC") != NULL ||
396 strstr(fname,".bmp" ) != NULL ||
397 strstr(fname,".BMP" ) != NULL ||
398 strstr(fname,".pbm" ) != NULL ||
399 strstr(fname,".PBM" ) != NULL ||
400 strstr(fname,".pgm" ) != NULL ||
401 strstr(fname,".PGM" ) != NULL ||
402 strstr(fname,".ppm" ) != NULL ||
403 strstr(fname,".PPM" ) != NULL ||
404 strstr(fname,".heic") != NULL ||
405 strstr(fname,".HEIC") != NULL ||
406 strstr(fname,".png" ) != NULL ||
407 strstr(fname,".PNG" ) != NULL ){
408
409 im = mri_read_stuff(fname) ; if( im != NULL ) RETURN(im) ;
410 }
411
412 /*-- 16 Aug 2006: AFNI dataset? --*/
413
414 if( strstr(fname,".HEAD") != NULL || strstr(fname,".BRIK") != NULL
415 || strstr(fname,".nii") != NULL ){
416 THD_3dim_dataset *dset = THD_open_dataset(fname) ;
417 if( dset != NULL ){
418 DSET_load(dset) ;
419 if( DSET_LOADED(dset) && DSET_datum_constant(dset) ){
420 im = mri_catvol_1D( dset->dblk->brick , 3 ) ;
421 if( im != NULL ){
422 im->dx = fabs(DSET_DX(dset)) ;
423 im->dy = fabs(DSET_DY(dset)) ;
424 im->dz = fabs(DSET_DZ(dset)) ;
425 }
426 }
427 DSET_delete(dset) ;
428 if( im != NULL ) RETURN(im) ;
429 }
430 }
431
432 /*-- check if file exists and is readable --*/
433
434 imfile = fopen( fname , "r" ) ;
435 if( imfile == NULL ){
436 ERROR_message("mri_read: couldn't open file %s" , fname ) ;
437 RETURN( NULL );
438 }
439
440 length = THD_filesize(fname) ; /* 22 Mar 2007 */
441
442 /*--- 03 Dec 2001: check for GEMS format file "IMGF" ---*/
443 /*[[[ Information herein from Medical Image Format FAQ ]]]*/
444
445 { char str[5]="AFNI" ;
446 int nx , ny , bpp , cflag , hdroff , extraskip=0 ;
447 rewind(imfile) ; fread(str,1,4,imfile) ; /* check for "IMGF" or "GEMS" */
448
449 if( str[0]=='G' && str[1]=='E' && str[2]=='M' && str[3]=='S' ){ /* 12 Feb 2004 */
450 char buf[4096]; int bb,cc; /* search for IMGF in 1st 4K */
451 rewind(imfile); cc=fread(buf,1,4096,imfile); cc-=4 ;
452 for( bb=4; bb < cc ; bb++ )
453 if( buf[bb]=='I' && buf[bb+1]=='M' && buf[bb+2]=='G' && buf[bb+3]=='F' ) break ;
454 if( bb < cc ){
455 fseek( imfile , (long)bb , SEEK_SET ) ; extraskip = bb ;
456 fread(str,1,4,imfile) ;
457 }
458 }
459
460 /* 12 Feb 2004: modified to allow for starting at extraskip */
461
462 if( str[0]=='I' && str[1]=='M' && str[2]=='G' && str[3]=='F' ){
463
464 fread( &skip , 4,1, imfile ) ; /* read next 5 ints */
465 fread( &nx , 4,1, imfile ) ;
466 fread( &ny , 4,1, imfile ) ;
467 fread( &bpp , 4,1, imfile ) ;
468 fread( &cflag, 4,1, imfile ) ;
469
470 if( nx < 0 || nx > 8192 ){ /* maybe have to byte swap 5 ints */
471 swap = 1 ; /* flag that we are swapping data */
472 swap_4(&skip); swap_4(&nx) ;
473 swap_4(&ny) ; swap_4(&bpp); swap_4(&cflag);
474 }
475 skip += extraskip ; /* location of image data in file */
476 if( nx < 0 || nx > 8192 || ny < 0 || ny > 8192 ) goto The_Old_Way ;
477 if( skip < 0 || skip >= length ) goto The_Old_Way ;
478 if( bpp != 16 || cflag != 1 ) goto The_Old_Way ;
479
480 /* make image space */
481
482 im = mri_new( nx , ny , MRI_short ) ;
483
484 /* try to read image auxiliary data as well (not mandatory) */
485
486 length = fseek( imfile , 148L+extraskip , SEEK_SET ) ; /* magic GEMS offset */
487 if( length == 0 ){
488 fread( &hdroff , 4,1 , imfile ) ; /* location of image header */
489 if( swap ) swap_4(&hdroff) ;
490 if( hdroff > 0 ){ /* read from image header */
491 float dx,dy,dz, dxx,dyy,dzz, xyz[9], zz ; int itr, ii,jj,kk, qq ;
492 static int nzoff=0 ;
493 static float zoff ;
494
495 /* get voxel grid sizes */
496
497 fseek( imfile , hdroff+26+extraskip , SEEK_SET ) ;
498 fread( &dzz , 4,1 , imfile ) ;
499
500 fseek( imfile , hdroff+50+extraskip , SEEK_SET ) ;
501 fread( &dxx , 4,1 , imfile ) ;
502 fread( &dyy , 4,1 , imfile ) ;
503
504 if( swap ){ swap_4(&dxx); swap_4(&dyy); swap_4(&dzz); }
505
506 /* save into image header [dw > 0 is signal that dx,dy,dz are OK] */
507
508 if( dxx > 0.01 && dyy > 0.01 && dzz > 0.01 ){
509 im->dx = dxx; im->dy = dyy; im->dz = dzz; im->dw = 1.0;
510 }
511
512 /* grid orientation: from 3 sets of LPI corner coordinates: */
513 /* xyz[0..2] = top left hand corner of image (TLHC) */
514 /* xyz[3..5] = top right hand corner of image (TRHC) */
515 /* xyz[6..8] = bottom right hand corner of image (BRHC) */
516 /* GEMS coordinate orientation here is LPI. */
517 /* Orientation is saved into global string MRILIB_orients. */
518 /* N.B.: AFNI coordinates are RAI orientation. */
519
520 fseek( imfile , hdroff+154+extraskip , SEEK_SET ) ;
521 fread( xyz , 4,9 , imfile ) ;
522 if( swap ) swap_fourbytes(9,xyz) ;
523
524 /* x-axis orientation */
525 /* ii determines which spatial direction is x-axis */
526 /* and is the direction that has the biggest change */
527 /* between the TLHC and TRHC */
528
529 dx = fabs(xyz[3]-xyz[0]) ; ii = 1 ;
530 dy = fabs(xyz[4]-xyz[1]) ; if( dy > dx ){ ii=2; dx=dy; }
531 dz = fabs(xyz[5]-xyz[2]) ; if( dz > dx ){ ii=3; }
532 dx = xyz[ii+2]-xyz[ii-1] ; if( dx < 0. ){ ii = -ii; }
533 switch( ii ){
534 case 1: MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; break;
535 case -1: MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; break;
536 case 2: MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; break;
537 case -2: MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; break;
538 case 3: MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; break;
539 case -3: MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; break;
540 default: MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; break;
541 }
542
543 /* y-axis orientation */
544 /* jj determines which spatial direction is y-axis */
545 /* and is the direction that has the biggest change */
546 /* between the BRHC and TRHC */
547
548 dx = fabs(xyz[6]-xyz[3]) ; jj = 1 ;
549 dy = fabs(xyz[7]-xyz[4]) ; if( dy > dx ){ jj=2; dx=dy; }
550 dz = fabs(xyz[8]-xyz[5]) ; if( dz > dx ){ jj=3; }
551 dx = xyz[jj+5]-xyz[jj+2] ; if( dx < 0. ){ jj = -jj; }
552 switch( jj ){
553 case 1: MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; break;
554 case -1: MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; break;
555 case 2: MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; break;
556 case -2: MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; break;
557 case 3: MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; break;
558 case -3: MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; break;
559 default: MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; break;
560 }
561
562 MRILIB_orients[6] = '\0' ; /* terminate orientation string */
563
564 kk = 6 - abs(ii)-abs(jj) ; /* which spatial direction is z-axis */
565 /* where 1=L-R, 2=P-A, 3=I-S */
566 zz = xyz[kk-1] ; /* z-coordinate of this slice (from TLHC) */
567
568 im->zo = zz ; /* 07 Aug 2002: save slice offset */
569
570 /* getting orientation of z-axis requires 2 images in a row -*/
571
572 if( nzoff == 0 ){ /* from 1st GEMS image */
573
574 zoff = zz ; /* save this for 2nd image calculation */
575 switch( kk ){ /* may be changed on second image */
576 case 1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
577 case 2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
578 case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
579 default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
580 }
581
582 } else if( nzoff == 1 ){ /* from 2nd GEMS image */
583
584 float qoff = zz - zoff ; /* vive la difference */
585 if( qoff < 0 ) kk = -kk ; /* kk determines z-axis orientation */
586
587 if( !use_MRILIB_slicespacing && qoff != 0.0 ){ /* 10 Jan 2004 */
588 use_MRILIB_slicespacing = 1 ;
589 MRILIB_slicespacing = fabs(qoff) ;
590 }
591
592 switch( kk ){
593 case 1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
594 case -1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
595 case 2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
596 case -2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
597 case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
598 case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
599 default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
600 }
601
602 /* save spatial offset of first slice */
603 /* [this needs to be positive in the direction of] */
604 /* [the -z axis, so may need to change its sign ] */
605
606 MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
607 if( kk == 1 || kk == 2 || kk == 3 ) MRILIB_zoff = -MRILIB_zoff ;
608
609 /* Same for x offset; [20 Dec 2001]
610 This must be at the middle of the TLHC voxel,
611 so we must move a little bit towards the TRHC edge;
612 We only use the result if the x-coordinate doesn't
613 change significantly between the TRHC and BRHC,
614 to avoid problems with oblique slices. */
615
616 qq = abs(ii) ;
617 MRILIB_xoff = ( xyz[qq-1]*(nx-0.5) + xyz[qq+2]*0.5 ) / nx ;
618 if( ii == 1 || ii == 2 || ii == 3 ) MRILIB_xoff = -MRILIB_xoff ;
619 use_MRILIB_xoff = ( fabs(xyz[qq+2]-xyz[qq+5]) < 0.01*dxx ) ;
620
621 /* Same for y offset;
622 This must be at the middle of the TLHC voxel,
623 so we must move a little bit towards the BRHC edge;
624 We only use the result if the y-coordinate doesn't
625 change significantly between the TLHC and TRHC. */
626
627 qq = abs(jj) ;
628 MRILIB_yoff = ( xyz[qq-1]*(ny-0.5) + xyz[qq+5]*0.5 ) / ny ;
629 if( jj == 1 || jj == 2 || jj == 3 ) MRILIB_yoff = -MRILIB_yoff ;
630 use_MRILIB_yoff = ( fabs(xyz[qq-1]-xyz[qq+2]) < 0.01*dyy ) ;
631 }
632 nzoff++ ; /* 3rd and later images don't count for z-orientation */
633
634 /* get TR, save into global variable */
635
636 if( MRILIB_tr <= 0.0 ){
637 fseek( imfile , hdroff+194+extraskip , SEEK_SET ) ;
638 fread( &itr , 4,1 , imfile ) ;
639 if( swap ) swap_4(&itr) ;
640 MRILIB_tr = im->dt = 1.0e-6 * itr ;
641 }
642 }
643 } /* end of trying to read image header */
644
645 goto Ready_To_Roll ; /* skip to the reading place */
646 }
647 } /* end of GEMS */
648
649 /*--- OK, do it the old way ---*/
650
651 The_Old_Way:
652
653 #if 0
654 MRILIB_orients[0] = '\0' ; MRILIB_zoff = MRILIB_tr = 0.0 ; /* 03 Dec 2001 */
655 #endif
656
657 switch( length ){
658
659 case 512: /* raw 16x16 short -- RWCox: 06 Dec 2001 */
660 im = mri_new( 16 , 16 , MRI_short ) ;
661 break ;
662
663 case 2048: /* raw 32x32 short -- RWCox: 19 Sep 2000 */
664 im = mri_new( 32 , 32 , MRI_short ) ;
665 break ;
666
667 case 4096: /* raw 64x64 byte -- RWC 3/21/95 */
668 im = mri_new( 64 , 64 , MRI_byte ) ;
669 break ;
670
671 case 8192: /* raw 64x64 short */
672 case 16096: /* with Signa 5.x header */
673 im = mri_new( 64 , 64 , MRI_short ) ;
674 skip = length - 8192 ;
675 break ;
676
677 #if 0
678 case 18432: /* raw 96x96 short */
679 im = mri_new( 96 , 96 , MRI_short ) ;
680 break ;
681 #endif
682
683 case 16384: /* raw 128x128 byte -- RWC 3/21/95 */
684 im = mri_new( 128 , 128 , MRI_byte ) ;
685 break ;
686
687 case 32768: /* raw 128x128 short */
688 case 40672: /* with Signa 5.x header */
689 im = mri_new( 128 , 128 , MRI_short ) ;
690 skip = length - 32768 ;
691 break ;
692
693 case 65536: /* raw 256x256 8-bit -- Matthew Belmonte March 1995 */
694 im = mri_new( 256 , 256 , MRI_byte ) ;
695 break ;
696
697 case 131072: /* raw 256x256 short */
698 case 138976: /* Signa 5.x */
699 case 145408: /* Signa 4.x */
700
701 im = mri_new( 256 , 256 , MRI_short ) ;
702 skip = length - 131072 ;
703 break ;
704
705 #if 0
706 case 262144: /* raw 256x256 float */
707 im = mri_new( 256 , 256 , MRI_float ) ;
708 break ;
709 #else
710 case 262144: /* raw 512x512 byte -- RWC 3/21/95 */
711 im = mri_new( 512 , 512 , MRI_byte ) ;
712 break ;
713
714 case 524288: /* raw 512x512 short -- RWC 3/21/95 */
715 im = mri_new( 512 , 512 , MRI_short ) ;
716 break ;
717
718 case 1048576: /* raw 1024x1024 byte -- RWC 3/21/95 */
719 im = mri_new( 1024 , 1024 , MRI_byte ) ;
720 break ;
721
722 case 2097152: /* raw 1024x1024 short -- RWC 3/21/95 */
723 im = mri_new( 1024 , 1024 , MRI_short ) ;
724 break ;
725 #endif
726
727 /** not a canonical length: try something else **/
728
729 default:
730 im = mri_try_mri( imfile , &skip ) ; /* Cox format */
731 if( im == NULL ) im = mri_try_7D ( imfile , &skip ) ; /* 7D format */
732 if( im == NULL ) im = mri_try_pgm( imfile , &skip ) ; /* PGM format */
733 if( im != NULL ) break ;
734
735 fclose( imfile ) ; /* close it, since we failed (so far) */
736
737 im = mri_read_ascii( fname ) ; /* list of ASCII numbers */
738 if( im != NULL ) RETURN( im );
739
740 im = mri_read_ppm( fname ) ; /* 15 Apr 1999 */
741 if( im != NULL ) RETURN( im );
742
743 im = mri_read_stuff( fname ) ; /* 22 Nov 2002 */
744 if( im != NULL ) RETURN( im );
745
746 ERROR_message("mri_read: do not recognize image file %s" , fname );
747 ERROR_message(" -- FWIW: length seen as %d bytes" , length ) ;
748 RETURN( NULL );
749 }
750
751 /*-- Actually read the data from disk --*/
752
753 Ready_To_Roll:
754
755 data = mri_data_pointer( im ) ;
756
757 length = fseek( imfile , skip , SEEK_SET ) ;
758 if( length != 0 ){
759 fprintf( stderr , "mri_read error in skipping in file %s\n" , fname ) ;
760 mri_free( im ) ;
761 RETURN( NULL );
762 }
763
764 length = fread( data , im->pixel_size , im->nvox , imfile ) ;
765 fclose( imfile ) ;
766
767 if( length != im->nvox ){
768 mri_free( im ) ;
769 fprintf( stderr , "couldn't read image data from file %s\n" , fname ) ;
770 RETURN( NULL );
771 }
772
773 mri_add_name( fname , im ) ;
774
775 /*-- 03 Dec 2001: maybe need to swap bytes --*/
776
777 if( swap ){
778 switch( im->pixel_size ){
779 default: break ;
780 case 2: swap_twobytes ( im->nvox, data ) ; break ; /* short */
781 case 4: swap_fourbytes( im->nvox, data ) ; break ; /* int, float */
782 case 8: swap_fourbytes( 2*im->nvox, data ) ; break ; /* complex */
783 }
784
785 im->was_swapped = 1 ; /* 07 Mar 2002 */
786 }
787
788 RETURN( im );
789 }
790
791
792 /******************************************************************/
793
794 /* GEMS 4.x image reading - rickr 03 Jun 2003 */
795
796 /*! Read a single 2D GEMS 4.x image.
797
798 \param filename is the name of the file to try to read
799 \return NULL if an image could not be read, otherwise
800 the address of a new MRI_IMAGE structure
801 */
802
mri_read_ge4(char * filename)803 MRI_IMAGE * mri_read_ge4( char * filename )
804 {
805 MRI_IMAGE * im;
806 ge4_header H;
807
808 ENTRY( "mri_read_ge4" );
809
810 if ( filename == NULL )
811 {
812 fprintf( stderr, "** mri_read_ge4 - missing filename\n" );
813 RETURN( NULL );
814 }
815
816 /* try to read image file - return with image */
817 if ( ge4_read_header( &H, filename, True ) != 0 )
818 RETURN( NULL );
819
820 /* these dimensions are fixed */
821 if ( (im = mri_new(256, 256, MRI_short)) == NULL )
822 {
823 free(H.image);
824 RETURN( NULL );
825 }
826
827 /* fill im struct with data from H */
828 im->zo = H.im_h.im_loc; /* this may well be incorrect */
829 im->dt = H.im_h.tr;
830 im->was_swapped = H.swap;
831
832 if ( ( H.ser_h.fov > 1.0 ) &&
833 ( H.ser_h.fov < 1000.0 ) &&
834 ( H.ser_h.scan_mat_x > 0 ) &&
835 ( H.ser_h.scan_mat_x < 1000 ) &&
836 ( H.ser_h.scan_mat_y > 0 ) &&
837 ( H.ser_h.scan_mat_y < 1000 ) )
838 {
839 /* attempt to set dx, dy and dz from these */
840
841 im->dx = 2 * H.ser_h.fov / H.ser_h.scan_mat_x;
842 im->dy = im->dx;
843 im->dz = 2 * H.ser_h.fov / H.ser_h.scan_mat_y;
844 im->dw = 1;
845 }
846
847 memcpy( mri_data_pointer(im), H.image, H.im_bytes );
848
849 mri_add_name( filename, im );
850
851 free(H.image); /* your services are no longer required */
852
853 RETURN( im );
854 }
855
856
857 /*********************************************************************/
858
859 #if 0
860 #define NUMSCAN(var) \
861 { while( isspace(ch) ) {ch = getc(imfile) ;} \
862 if(ch == '#') do{ch = getc(imfile) ;}while(ch != '\n' && ch != EOF) ; \
863 if(ch =='\n') ch = getc(imfile) ; \
864 for( nch=0 ; isdigit(ch) ; nch++,ch=getc(imfile) ) {buf[nch] = ch ;} \
865 buf[nch]='\0'; \
866 var = strtol( buf , NULL , 10 ) ; }
867 #else
868
869 /*! Skip comments in a PPM file. */
870
871 #define SKIPCOM \
872 {if(ch == '#') do{ch = getc(imfile) ;}while(ch != '\n' && ch != EOF);}
873
874 /*! Scan for a number in the imfile stream. */
875
876 #define NUMSCAN(var) \
877 { SKIPCOM ; \
878 while( ch!=EOF && !isdigit(ch) ){ch = getc(imfile); SKIPCOM; } \
879 for( nch=0 ; isdigit(ch) ; nch++,ch=getc(imfile) ) {buf[nch] = ch ;} \
880 buf[nch]='\0'; \
881 var = strtol( buf , NULL , 10 ) ; }
882 #endif
883
884 /*! Try to read an file in the "Cox MRI" format.
885
886 \param imfile is a pointer to an open FILE.
887 \param skip is a pointer to an int that will be set to the number
888 of bytes to skip from the file start to find the image data
889 \return NULL if the file doesn't work for "Cox MRI" format;
890 otherwise, the return is a pointer to an MRI_IMAGE ready
891 to have its data read from imfile.
892 */
893
mri_try_mri(FILE * imfile,int * skip)894 static MRI_IMAGE *mri_try_mri( FILE *imfile , int *skip )
895 {
896 int ch , nch , nx,ny,imcode ;
897 char buf[64] ;
898 MRI_IMAGE *im ;
899
900 ENTRY("mri_try_mri") ;
901
902 fseek( imfile , 0 , SEEK_SET ) ; /* rewind file */
903
904 ch = getc( imfile ) ; if( ch != 'M' ) RETURN( NULL ); /* check for MRI */
905 ch = getc( imfile ) ; if( ch != 'R' ) RETURN( NULL );
906 ch = getc( imfile ) ; if( ch != 'I' ) RETURN( NULL );
907
908 /* magic MRI found, so read numbers */
909
910 ch = getc(imfile) ;
911
912 NUMSCAN(imcode) ;
913 NUMSCAN(nx) ; if( nx <= 0 ) RETURN( NULL );
914 NUMSCAN(ny) ; if( ny <= 0 ) RETURN( NULL );
915
916 *skip = ftell(imfile) ;
917 im = mri_new( nx , ny , imcode ) ;
918 RETURN( im );
919 }
920
921 /**************************************************************************
922 7D format: MRn kind n-dimensions data, where 'n' = 1-7.
923 ***************************************************************************/
924
925 /*! Try to read a "Cox nD MRI" image file (fat chance). */
926
mri_try_7D(FILE * imfile,int * skip)927 static MRI_IMAGE *mri_try_7D( FILE *imfile , int *skip )
928 {
929 int ch , nch , nx,ny,nz,nt,nu,nv,nw , imcode , ndim ;
930 char buf[64] ;
931 MRI_IMAGE *im ;
932
933 ENTRY("mri_try_7D") ;
934
935 fseek( imfile , 0 , SEEK_SET ) ; /* rewind file */
936
937 ch = getc( imfile ) ; if( ch != 'M' ) RETURN( NULL ); /* check for MR[1-7] */
938 ch = getc( imfile ) ; if( ch != 'R' ) RETURN( NULL );
939 ch = getc( imfile ) ;
940 switch( ch ){
941 default: RETURN( NULL ); /* not what I expected */
942
943 case '1': ndim = 1 ; break ;
944 case '2': ndim = 2 ; break ;
945 case '3': ndim = 3 ; break ;
946 case '4': ndim = 4 ; break ;
947 case '5': ndim = 5 ; break ;
948 case '6': ndim = 6 ; break ;
949 case '7': ndim = 7 ; break ;
950 }
951 /* magic MR? found, so read numbers */
952
953 ch = getc(imfile) ;
954 NUMSCAN(imcode) ;
955
956 nx = ny = nz = nt = nu = nv = nw = 1 ;
957
958 NUMSCAN(nx) ; if( nx <= 0 ) RETURN( NULL );
959 if( ndim > 1 ){ NUMSCAN(ny) ; if( ny <= 0 ) RETURN( NULL ); }
960 if( ndim > 2 ){ NUMSCAN(nz) ; if( nz <= 0 ) RETURN( NULL ); }
961 if( ndim > 3 ){ NUMSCAN(nt) ; if( nt <= 0 ) RETURN( NULL ); }
962 if( ndim > 4 ){ NUMSCAN(nu) ; if( nu <= 0 ) RETURN( NULL ); }
963 if( ndim > 5 ){ NUMSCAN(nv) ; if( nv <= 0 ) RETURN( NULL ); }
964 if( ndim > 6 ){ NUMSCAN(nw) ; if( nw <= 0 ) RETURN( NULL ); }
965
966 *skip = ftell(imfile) ;
967 im = mri_new_7D_generic( nx,ny,nz,nt,nu,nv,nw , imcode , TRUE ) ;
968 RETURN( im );
969 }
970
971 /*********************************************************************/
972
973 /*! Try to read a raw PGM format image file.
974
975 \param imfile is a pointer to an open FILE
976 \param skip is a pointer to an int; *skip will be set to the
977 byte offset at which to start reading data
978 \return A pointer to an MRI_IMAGE ready to have its data read in
979 (if the file is a PGM file), or NULL.
980 */
981
mri_try_pgm(FILE * imfile,int * skip)982 static MRI_IMAGE *mri_try_pgm( FILE *imfile , int *skip )
983 {
984 int ch , nch , nx,ny,maxval ;
985 char buf[64] ;
986 MRI_IMAGE *im ;
987
988 ENTRY("mri_try_pgm") ;
989
990 fseek( imfile , 0 , SEEK_SET ) ; /* rewind file */
991
992 ch = getc( imfile ) ; if( ch != 'P' ) RETURN(NULL); /* check for magic */
993 ch = getc( imfile ) ; if( ch != '5' ) RETURN(NULL);
994
995 /* magic P5 found, so read numbers */
996
997 ch = getc(imfile) ;
998
999 NUMSCAN(nx) ; if( nx <= 0 ) RETURN(NULL);
1000 NUMSCAN(ny) ; if( ny <= 0 ) RETURN(NULL);
1001 NUMSCAN(maxval) ; if( maxval <= 0 || maxval > 255 ) RETURN(NULL);
1002
1003 *skip = ftell(imfile) ;
1004 im = mri_new( nx , ny , MRI_byte ) ;
1005 RETURN(im);
1006 }
1007
1008 /*--------------------------------------------------------------
1009 Read a pile of images from one file.
1010 Modified 4/4/95 to read short or byte data.
1011 Modified 10/02/95 to allow byte swapping with 3Ds:
1012 Modified 11/06/95 to allow float images with 3Df:
1013 and to allow int images with 3Di:
1014 and to allow complex images with 3Dc:
1015 Modified 16 Apr 2002 to allow RGB input with 3Dr:
1016
1017 [N.B.: if this routine is altered, don't forget mri_imcount!]
1018 ----------------------------------------------------------------*/
1019
1020 /*! Read one or more 2D slices from a "3D:" formatted image file. */
1021
mri_read_3D(char * tname)1022 MRI_IMARR * mri_read_3D( char *tname )
1023 {
1024 int hglobal , himage , nx , ny , nz ;
1025 char fname[256] , buf[512] ;
1026 int ngood , kim , datum_type , datum_len , swap ;
1027 MRI_IMARR *newar ;
1028 MRI_IMAGE *newim ;
1029 void *imar ;
1030 FILE *imfile ;
1031 long long length , nneed , koff , hglob ; /* 22 Mar 2007 */
1032
1033 ENTRY("mri_read_3D") ;
1034
1035 /*** get info from 3D tname ***/
1036
1037 if( tname == NULL || strlen(tname) < 10 ) RETURN(NULL) ;
1038
1039 switch( tname[2] ){ /* allow for 3D: or 3Ds: or 3Db:, etc */
1040
1041 default:
1042 case ':':
1043 ngood = sscanf( tname , "3D:%d:%d:%d:%d:%d:%s" ,
1044 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1045
1046 swap = 0 ;
1047 datum_type = MRI_short ;
1048 datum_len = sizeof(short) ; /* better be 2 */
1049 break ;
1050
1051 case 's':
1052 ngood = sscanf( tname , "3Ds:%d:%d:%d:%d:%d:%s" ,
1053 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1054
1055 swap = 1 ;
1056 datum_type = MRI_short ;
1057 datum_len = sizeof(short) ; /* better be 2 */
1058 break ;
1059
1060 case 'b':
1061 ngood = sscanf( tname , "3Db:%d:%d:%d:%d:%d:%s" ,
1062 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1063
1064 swap = 0 ;
1065 datum_type = MRI_byte ;
1066 datum_len = sizeof(byte) ; /* better be 1 */
1067 break ;
1068
1069 case 'f':
1070 ngood = sscanf( tname , "3Df:%d:%d:%d:%d:%d:%s" ,
1071 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1072
1073 swap = 0 ;
1074 datum_type = MRI_float ;
1075 datum_len = sizeof(float) ; /* better be 4 */
1076 break ;
1077
1078 case 'd': /* 06 Feb 2003 */
1079 ngood = sscanf( tname , "3Dd:%d:%d:%d:%d:%d:%s" ,
1080 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1081
1082 swap = 0 ;
1083 datum_type = MRI_double ;
1084 datum_len = sizeof(double) ; /* better be 8 */
1085 break ;
1086
1087 case 'i':
1088 ngood = sscanf( tname , "3Di:%d:%d:%d:%d:%d:%s" ,
1089 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1090
1091 swap = 0 ;
1092 datum_type = MRI_int ;
1093 datum_len = sizeof(int) ; /* better be 4 */
1094 break ;
1095
1096 case 'c':
1097 ngood = sscanf( tname , "3Dc:%d:%d:%d:%d:%d:%s" ,
1098 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1099
1100 swap = 0 ;
1101 datum_type = MRI_complex ;
1102 datum_len = sizeof(complex) ; /* better be 8 */
1103 break ;
1104
1105 case 'r':
1106 ngood = sscanf( tname , "3Dr:%d:%d:%d:%d:%d:%s" ,
1107 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1108
1109 swap = 0 ;
1110 datum_type = MRI_rgb ;
1111 datum_len = 3*sizeof(byte) ; /* better be 3 */
1112 break ;
1113 }
1114
1115 if( ngood < 6 || himage < 0 ||
1116 nx <= 0 || ny <= 0 || nz <= 0 ||
1117 strlen(fname) <= 0 ) RETURN(NULL); /* bad info */
1118
1119 /*** 06 Mar 2001: special case of fname ***/
1120
1121 if( strcmp(fname,"ALLZERO") == 0 ){
1122 INIT_IMARR(newar) ;
1123 for( kim=0 ; kim < nz ; kim++ ){
1124 newim = mri_new( nx , ny , datum_type ) ;
1125 imar = mri_data_pointer( newim ) ;
1126 memset( imar , 0 , newim->nvox * newim->pixel_size ) ;
1127 sprintf( buf , "%s#%d" , fname,kim ) ;
1128 mri_add_name( buf , newim ) ;
1129 ADDTO_IMARR(newar,newim) ;
1130 }
1131 RETURN(newar);
1132 }
1133
1134 /*** open the input file and position it ***/
1135
1136 imfile = fopen( fname , "r" ) ;
1137 if( imfile == NULL ){
1138 ERROR_message("mri_read_3D: couldn't open file %s" , fname ) ;
1139 RETURN(NULL);
1140 }
1141
1142 length = THD_filesize(fname) ; /* 22 Mar 2007 */
1143
1144 /** 13 Apr 1999: modified to allow actual hglobal < -1
1145 as long as hglobal+himage >= 0 **/
1146
1147 hglob = hglobal ;
1148 if( hglob == -1 || hglob+himage < 0 ){
1149 hglob = length - (datum_len*nx*ny+himage) * (long long)nz ;
1150 if( hglob < 0 ) hglob = 0 ;
1151 }
1152
1153 nneed = hglob + (datum_len*nx*ny+himage) * (long long)nz ;
1154 if( length < nneed ){
1155 ERROR_message(
1156 "mri_read_3D: file %s is %lld bytes long but must be at least %lld bytes long\n"
1157 " for hglobal=%lld himage=%d nx=%d ny=%d nz=%d and voxel=%d bytes\n",
1158 fname,length,nneed,hglob,himage,nx,ny,nz,datum_len ) ;
1159 fclose( imfile ) ;
1160 RETURN(NULL);
1161 }
1162
1163 /*** read images from the file ***/
1164
1165 INIT_IMARR(newar) ;
1166
1167 for( kim=0 ; kim < nz ; kim++ ){
1168 koff = hglob + (kim+1)*himage + datum_len*nx*ny * (long long)kim ;
1169 fseeko( imfile, (off_t)koff, SEEK_SET ) ; /* 22 Mar 2007: fseek->fseeko */
1170
1171 newim = mri_new( nx , ny , datum_type ) ;
1172 imar = mri_data_pointer( newim ) ;
1173 (void)fread( imar , datum_len , nx * ny , imfile ) ;
1174 if( swap ){
1175 mri_swapbytes( newim ) ;
1176 newim->was_swapped = 1 ; /* 07 Mar 2002 */
1177 }
1178
1179 if( nz == 1 ) mri_add_name( fname , newim ) ;
1180 else {
1181 sprintf( buf , "%s#%d" , fname,kim ) ;
1182 mri_add_name( buf , newim ) ;
1183 }
1184
1185 ADDTO_IMARR(newar,newim) ;
1186 }
1187
1188 fclose(imfile) ;
1189 RETURN(newar);
1190 }
1191
1192 /*------------------------------------------------------------------------------*/
1193
1194 /*! Read one or more 2D images from a file.
1195
1196 This function is the main point of input for to3d.c.
1197 \param fname is the name of the file to read. This file
1198 might be in one of these formats:
1199 - "3D:" format (implicitly or explicitly)
1200 - "3A:" format
1201 - *.hdr (ANALYZE 2D-4D) format
1202 - *.ima (Siemens 2D array) format
1203 - I.* (GEMS) format
1204 - PGM format
1205 - PPM format
1206 - GIF, TIFF, JPEG, BMP, PNG formats (thru filters)
1207 - List of ASCII numbers
1208 - pre-defined 2D file size in mri_read()
1209 - "Cox MRI" (if this is what you have, God help you, no one else can)
1210
1211 \return A pointer to an array of 2D images. If nothing
1212 could be read, NULL is returned.
1213 */
1214
mri_read_file(char * fname)1215 MRI_IMARR * mri_read_file( char * fname )
1216 {
1217 MRI_IMARR *newar=NULL ;
1218 MRI_IMAGE *newim ;
1219 char *new_fname ;
1220 int tried_dicom=0 ;
1221
1222 ENTRY("mri_read_file") ;
1223
1224 /* convert fname to new_fname, based on environment */
1225
1226 new_fname = imsized_fname( fname ) ;
1227 if( new_fname == NULL ) RETURN( NULL );
1228
1229 /* input method is based on filename */
1230
1231 if( strlen(new_fname) > 9 &&
1232 new_fname[0] == '3' &&
1233 new_fname[1] == 'D' &&
1234 (new_fname[2] == ':' || new_fname[3] == ':') ){
1235
1236 newar = mri_read_3D( new_fname ) ; /* read from a 3D: file */
1237
1238 } else if( strlen(new_fname) > 9 &&
1239 new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
1240
1241 newar = mri_read_3A( new_fname ) ; /* from a 3A: file */
1242
1243 /*-- from a dataset? [10 Dec 2007] --*/
1244
1245 } else if( strstr(fname,".HEAD") != NULL || strstr(fname,".nii") != NULL ){
1246
1247 THD_3dim_dataset *dset = THD_open_dataset(fname) ;
1248 if( dset != NULL ){
1249 int ii,jj ; MRI_IMAGE *qim ; void *qar ; MRI_IMARR *qimar ;
1250 DSET_load(dset) ; INIT_IMARR(newar) ;
1251 for( ii=0 ; ii < DSET_NVALS(dset) ; ii++ ){
1252 qim = DSET_BRICK(dset,ii) ; qar = mri_data_pointer(qim) ;
1253 if( qim != NULL && qar != NULL ){
1254 qimar = mri_to_imarr( qim ) ;
1255 if( qimar != NULL ){
1256 for( jj=0 ; jj < IMARR_COUNT(qimar) ; jj++ )
1257 ADDTO_IMARR(newar,IMARR_SUBIM(qimar,jj)) ;
1258 FREE_IMARR(qimar) ;
1259 }
1260 }
1261 }
1262 DSET_delete(dset) ; RETURN(newar) ;
1263 }
1264
1265 } else if( check_dicom_magic_num( new_fname ) ) { /* 10 Aug 2004 */
1266
1267 newar = mri_read_dicom( new_fname ); tried_dicom=2 ;
1268
1269 } else if( strstr(new_fname,".hdr") != NULL ||
1270 strstr(new_fname,".HDR") != NULL ){ /* 05 Feb 2001 */
1271
1272 newar = mri_read_analyze75( new_fname ) ; /* ANALYZE .hdr/.img filepair */
1273
1274 } else if( strstr(new_fname,".ima") != NULL ||
1275 strstr(new_fname,".IMA") != NULL ){ /* 12 Mar 2001 */
1276
1277 newar = mri_read_siemens( new_fname ) ; /* Siemens file */
1278
1279 } else if( strncmp(new_fname,"I.",2) == 0 || /* GE I.* files */
1280 strstr(new_fname,"/I.") != NULL ||
1281 strstr(new_fname,".ppm") != NULL || /* raw PPM or PGM files */
1282 strstr(new_fname,".pgm") != NULL ||
1283 strstr(new_fname,".pnm") != NULL ||
1284 strstr(new_fname,".PPM") != NULL ||
1285 strstr(new_fname,".PNM") != NULL ||
1286 strstr(new_fname,".PGM") != NULL ){ /* 05 Nov 2002 */
1287
1288 newim = mri_read( new_fname ) ; /* read from a 2D file with 1 slice */
1289
1290 if ( newim == NULL ) /* GEMS 4.x - 03 Jun 2003 [rickr] */
1291 newim = mri_read_ge4( new_fname ) ;
1292
1293 if( newim != NULL ){
1294 INIT_IMARR(newar) ;
1295 ADDTO_IMARR(newar,newim) ;
1296 }
1297
1298 } else if( strncmp(new_fname,"i.",2) == 0 || /* GEMS 4.x i.* files */
1299 strstr(new_fname,"/i.") != NULL ){ /* 03 Jun 2003 [rickr] */
1300
1301 newim = mri_read_ge4( new_fname ) ; /* 2D file with 1 slice */
1302
1303 if( newim != NULL ){
1304 INIT_IMARR(newar) ;
1305 ADDTO_IMARR(newar,newim) ;
1306 }
1307
1308 } else if( strstr(new_fname,".jpg" ) != NULL || /* various formats */
1309 strstr(new_fname,".JPG" ) != NULL || /* that we convert */
1310 strstr(new_fname,".jpeg") != NULL || /* to PPG/PGM using */
1311 strstr(new_fname,".JPEG") != NULL || /* external filters */
1312 strstr(new_fname,".gif" ) != NULL ||
1313 strstr(new_fname,".GIF" ) != NULL ||
1314 strstr(new_fname,".tif" ) != NULL ||
1315 strstr(new_fname,".TIF" ) != NULL ||
1316 strstr(new_fname,".tiff") != NULL ||
1317 strstr(new_fname,".TIFF") != NULL ||
1318 strstr(new_fname,".heic") != NULL ||
1319 strstr(new_fname,".HEIC") != NULL ||
1320 strstr(new_fname,".bmp" ) != NULL ||
1321 strstr(new_fname,".BMP" ) != NULL ||
1322 strstr(new_fname,".pbm" ) != NULL ||
1323 strstr(new_fname,".PBM" ) != NULL ||
1324 strstr(new_fname,".png" ) != NULL ||
1325 strstr(new_fname,".PNG" ) != NULL ){ /* 22 Nov 2002 */
1326
1327 newim = mri_read_stuff( new_fname ) ;
1328 if( newim != NULL ){
1329 INIT_IMARR(newar) ;
1330 ADDTO_IMARR(newar,newim) ;
1331 }
1332
1333 } else if( strstr(new_fname,".mpg" ) != NULL || /* 03 Dec 2003 */
1334 strstr(new_fname,".MPG" ) != NULL || /* read MPEGs */
1335 strstr(new_fname,".mpeg") != NULL ||
1336 strstr(new_fname,".MPEG") != NULL ){
1337
1338 newar = mri_read_mpeg( new_fname ) ; /* cf. mri_read_mpeg.c */
1339 }
1340
1341 /** failed to read anything? try DICOM format (doesn't have a fixed suffix) **/
1342 /* 05 May 2003 added option to try DICOM last KRH */
1343
1344 if( newar == NULL ){
1345
1346 if ( !AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
1347 if( !tried_dicom ){
1348 newar = mri_read_dicom( new_fname ) ; tried_dicom = 1 ;
1349 }
1350 }
1351
1352 /** if DICOM failed, try a 2D slice file, hope for the best **/
1353
1354 if( newar == NULL && tried_dicom != 2 ){
1355 newim = mri_read( new_fname ) ;
1356 if( newim == NULL ){ free(new_fname); RETURN( NULL ); } /* give up */
1357 INIT_IMARR(newar) ;
1358 ADDTO_IMARR(newar,newim) ;
1359 }
1360
1361 if( newar == NULL && AFNI_yesenv("AFNI_TRY_DICOM_LAST") ){
1362 if( !tried_dicom ){
1363 newar = mri_read_dicom( new_fname ) ; tried_dicom = 1 ;
1364 }
1365 }
1366 }
1367
1368 free(new_fname) ; /* done with the mangled filename */
1369
1370 /* 07 Mar 2002: add fname to the images, if needed */
1371
1372 if( newar != NULL && newar->num > 0 ){
1373 int ii ;
1374 for( ii=0 ; ii < newar->num ; ii++ ){
1375 newim = IMARR_SUBIM(newar,ii) ;
1376 if( newim != NULL && newim->fname == NULL )
1377 newim->fname = strdup(fname) ;
1378 }
1379 }
1380
1381 RETURN( newar );
1382 }
1383
1384 /*-----------------------------------------------------------------*/
1385
1386 /*! Like mri_read_file(), but will only return 1 2D image.
1387
1388 If the input file has more than 1 slice, or cannot be read,
1389 then NULL is returned.
1390 */
1391
mri_read_just_one(char * fname)1392 MRI_IMAGE * mri_read_just_one( char * fname )
1393 {
1394 MRI_IMARR * imar ;
1395 MRI_IMAGE * im ;
1396 char * new_fname ;
1397
1398 ENTRY("mri_read_just_one") ;
1399
1400 new_fname = imsized_fname( fname ) ;
1401 if( new_fname == NULL ) RETURN( NULL );
1402
1403 imar = mri_read_file( new_fname ) ; free(new_fname) ;
1404 if( imar == NULL ) RETURN( NULL );
1405 if( imar->num != 1 ){ DESTROY_IMARR(imar) ; RETURN( NULL ); }
1406 im = IMAGE_IN_IMARR(imar,0) ;
1407 FREE_IMARR(imar) ;
1408 RETURN( im );
1409 }
1410
1411 /*-----------------------------------------------------------------
1412 return a count of how many 2D images will be read from this file
1413 -------------------------------------------------------------------*/
1414
1415 /*! Return a count of how many 2D images are in a file.
1416
1417 Used by to3d.c to figure out how many slices will be read
1418 later using mri_read_file(). Return value is 0 if the images
1419 can't be counted. If you add a new file type to mri_read_file(),
1420 then you need to modify this function as well!
1421 */
1422
1423 static int mri_imcount_analyze75( char * ) ; /* prototype */
1424 static int mri_imcount_siemens( char * ) ;
1425
mri_imcount(char * tname)1426 int mri_imcount( char *tname )
1427 {
1428 int hglobal , himage , nx , ny , nz , ngood ;
1429 char fname[256]="\0" ;
1430 char *new_fname ;
1431
1432 ENTRY("mri_imcount") ;
1433
1434 if( tname == NULL ) RETURN( 0 );
1435 new_fname = imsized_fname( tname ) ;
1436 if( new_fname == NULL ) RETURN( 0 );
1437
1438 /*** a 3D dataset [06 Jan 2011] ***/
1439
1440 if( strstr(new_fname,".HEAD") != NULL || strstr(new_fname,".nii") != NULL ){
1441 THD_3dim_dataset *dset = THD_open_dataset(new_fname) ;
1442 if( dset != NULL ){
1443 nz = DSET_NZ(dset) * DSET_NVALS(dset) ; DSET_delete(dset) ;
1444 RETURN(nz) ;
1445 }
1446 }
1447
1448 /*** a 3D filename ***/
1449
1450 if( strlen(new_fname) > 9 && new_fname[0] == '3' && new_fname[1] == 'D' &&
1451 (new_fname[2] == ':' || new_fname[3] == ':') ){
1452 /* check for ':', too 3 Jan 2005 [rickr] */
1453 switch( new_fname[2] ){
1454
1455 default:
1456 case ':':
1457 ngood = sscanf( new_fname , "3D:%d:%d:%d:%d:%d:%s" ,
1458 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1459 break ;
1460
1461 case 's':
1462 ngood = sscanf( new_fname , "3Ds:%d:%d:%d:%d:%d:%s" ,
1463 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1464 break ;
1465
1466 case 'b':
1467 ngood = sscanf( new_fname , "3Db:%d:%d:%d:%d:%d:%s" ,
1468 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1469 break ;
1470
1471 case 'f':
1472 ngood = sscanf( new_fname , "3Df:%d:%d:%d:%d:%d:%s" ,
1473 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1474 break ;
1475
1476 case 'd': /* 06 Feb 2003 */
1477 ngood = sscanf( new_fname , "3Dd:%d:%d:%d:%d:%d:%s" ,
1478 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1479 break ;
1480
1481 case 'i':
1482 ngood = sscanf( new_fname , "3Di:%d:%d:%d:%d:%d:%s" ,
1483 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1484 break ;
1485
1486 case 'c':
1487 ngood = sscanf( new_fname , "3Dc:%d:%d:%d:%d:%d:%s" ,
1488 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1489 break ;
1490
1491 case 'r':
1492 ngood = sscanf( new_fname , "3Dr:%d:%d:%d:%d:%d:%s" ,
1493 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
1494 break ;
1495 }
1496
1497 free( new_fname ) ;
1498 if( ngood < 6 || himage < 0 ||
1499 nx <= 0 || ny <= 0 || nz <= 0 ||
1500 strlen(fname) <= 0 ) RETURN( 0 );
1501 else RETURN( nz );
1502 }
1503
1504 /*** a 3A filename ***/
1505
1506 if( strlen(new_fname) > 9 &&
1507 new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
1508
1509 switch( new_fname[2] ){
1510
1511 default: ngood = 0 ; break ;
1512
1513 case 's':
1514 ngood = sscanf( new_fname, "3As:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
1515 break ;
1516
1517 case 'b':
1518 ngood = sscanf( new_fname, "3Ab:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
1519 break ;
1520
1521 case 'f':
1522 ngood = sscanf( new_fname, "3Af:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
1523 break ;
1524 }
1525
1526 free( new_fname ) ;
1527 if( ngood < 4 || nx <= 0 || ny <= 0 || nz <= 0 || strlen(fname) <= 0 ) RETURN( 0 );
1528 else RETURN( nz );
1529 }
1530
1531 /*** 05 Feb 2001: deal with ANALYZE .hdr files ***/
1532
1533 if( strstr(new_fname,".hdr") != NULL ||
1534 strstr(new_fname,".HDR") != NULL ){
1535
1536 nz = mri_imcount_analyze75( new_fname ) ;
1537 if( nz > 0 ){ free(new_fname); RETURN(nz); }
1538 }
1539
1540 if( strstr(new_fname,".ima") != NULL ||
1541 strstr(new_fname,".IMA") != NULL ){ /* 12 Mar 2001 */
1542
1543 nz = mri_imcount_siemens( new_fname ) ;
1544 if( nz > 0 ){ free(new_fname); RETURN(nz); }
1545 }
1546
1547 if( strstr(new_fname,".mpg" ) != NULL || /* 03 Dec 2003 */
1548 strstr(new_fname,".MPG" ) != NULL ||
1549 strstr(new_fname,".mpeg") != NULL ||
1550 strstr(new_fname,".MPEG") != NULL ){
1551
1552 nz = mri_imcount_mpeg( new_fname ) ;
1553 if( nz > 0 ){ free(new_fname); RETURN(nz); }
1554 }
1555
1556 /*** 19 Jul 2002: see if it is a DICOM file ***/
1557
1558 mri_dicom_seterr(0) ;
1559 nz = mri_imcount_dicom( new_fname ) ; /* cf. mri_read_dicom.c */
1560 mri_dicom_seterr(1) ;
1561 if( nz > 0 ){ free(new_fname); RETURN(nz); }
1562
1563 /*** not recognized ***/
1564
1565 free(new_fname) ; RETURN(1) ; /* assume it has 1 image in it, somewhere */
1566 }
1567
1568 /*--------------------------------------------------------------*/
1569
1570 /*! Like mri_read_file(), but returns images from many files.
1571
1572 \param nf = Number of file names
1573 \param fn = Array of file name strings
1574 \return An array of 2D images (NULL if nothing was found)
1575
1576 Added 07 Mar 1995
1577 */
1578
mri_read_many_files(int nf,char * fn[])1579 MRI_IMARR * mri_read_many_files( int nf , char * fn[] )
1580 {
1581 MRI_IMARR * newar , * outar ;
1582 int kf , ii ;
1583
1584 ENTRY("mri_read_many_files") ;
1585
1586 if( nf <= 0 ) RETURN( NULL ); /* no inputs! */
1587 INIT_IMARR(outar) ; /* initialize output array */
1588
1589 for( kf=0 ; kf < nf ; kf++ ){
1590 newar = mri_read_file( fn[kf] ) ; /* read all images in this file */
1591
1592 if( newar == NULL ){ /* none? flush the output array! */
1593 fprintf(stderr,"cannot read images from file %s\n",fn[kf]) ;
1594 for( ii=0 ; ii < outar->num ; ii++ ) mri_free(outar->imarr[ii]) ;
1595 FREE_IMARR(outar) ;
1596 RETURN( NULL );
1597 }
1598
1599 for( ii=0 ; ii < newar->num ; ii++ ) /* move images to output array */
1600 ADDTO_IMARR( outar , newar->imarr[ii] ) ;
1601
1602 FREE_IMARR(newar) ; /* don't need this no more */
1603 }
1604 RETURN( outar );
1605 }
1606
1607 /*! Like mri_read_many_files(), but forces images to a certain resolution.
1608
1609 \param nf = Number of file names
1610 \param fn = Array of file name strings
1611 \param nx = number of pixels
1612 \param ny in x and y directions
1613 if nx is negative, then nx and ny are set
1614 to be the dimensions of the very first image
1615 read.
1616 if ny is negative then resample but do not modify
1617 the aspect ratio of the image. After resampling,
1618 padding is used to achieve final pixel count.
1619 \param pval padding value, if padding is to be used.
1620 \return An array of 2D images (NULL if nothing was found)
1621
1622 Added Jan 07
1623 */
mri_read_resamp_many_files(int nf,char * fn[],int nxnew,int nynew,byte pval)1624 MRI_IMARR * mri_read_resamp_many_files( int nf, char * fn[] , int nxnew,
1625 int nynew, byte pval)
1626 {
1627 MRI_IMARR * newar , * outar ;
1628 int kf , ii, nxi, nyi , keepaspect , bot;
1629 MRI_IMAGE * bim, *qim, *imin, *zim;
1630
1631 ENTRY("mri_read_resamp_many_files") ;
1632
1633 if( nf <= 0 ) RETURN( NULL ); /* no inputs! */
1634 INIT_IMARR(outar) ; /* initialize output array */
1635
1636 if (nynew < 0) {
1637 keepaspect = 1;
1638 nynew = -nynew;
1639 } else {
1640 keepaspect = 0;
1641 }
1642
1643 for( kf=0 ; kf < nf ; kf++ ){
1644 newar = mri_read_file( fn[kf] ) ; /* read all images in this file */
1645
1646 if( newar == NULL ){ /* none? flush the output array! */
1647 fprintf(stderr,"cannot read images from file %s\n",fn[kf]) ;
1648 for( ii=0 ; ii < outar->num ; ii++ ) mri_free(outar->imarr[ii]) ;
1649 FREE_IMARR(outar) ;
1650 RETURN( NULL );
1651 }
1652 if (nxnew < 0 && kf == 0) { /* set dimensions based on 1st image */
1653 nxnew = newar->imarr[0]->nx;
1654 nynew = newar->imarr[0]->ny;
1655 }
1656
1657 for( ii=0 ; ii < newar->num ; ii++ ) {/* move images to output array */
1658 imin = newar->imarr[ii];
1659 nxi = imin->nx;
1660 nyi = imin->ny;
1661 if (nxi != nxnew || nyi != nynew) { /* resampling needed (adapted from galler.c)*/
1662 float fx , fy , fxm ;
1663 fx = nxnew / (float)nxi ; fy = nynew / (float)nyi ;
1664 fxm = MIN(fx,fy) ;
1665 /* fprintf(stderr,"Resizing from %dx%d to %dx%d.\n fx = %.3f\n", nxi, nyi, nxnew, nynew, fxm); */
1666 if( fxm < 0.95f ){
1667 float sigma = 0.3456789f/fxm ;
1668 /* fprintf(stderr,"sigma %f\n", sigma); */
1669 if (imin->kind == MRI_rgb) {
1670 bim = mri_rgb_blur2D( sigma , imin ) ;
1671 } else {
1672 bim = mri_byte_blur2D( sigma , imin ) ;
1673 }
1674 } else bim = imin ;
1675 if (keepaspect && fx != fy) {
1676 if (fx < fy) {
1677 qim = mri_resize(bim, nxnew, (int)(fx*nyi));
1678 /* fprintf(stderr,"qim X now %dx%d\n", qim->nx, qim->ny); */
1679 bot = (nynew - (int)(fx*nyi))/2;
1680 zim = mri_valpad_2D( 0 , 0 ,
1681 bot, nynew-(int)(fx*nyi)-bot, qim, pval);
1682 if (qim != zim) mri_free(qim) ;
1683 qim = zim; zim = NULL;
1684 /* fprintf(stderr,"qim X padded %dx%d, bot=%d\n",
1685 qim->nx, qim->ny, bot); */
1686 } else {
1687 qim = mri_resize(bim, (int)(fy*nxi), nynew);
1688 /* fprintf(stderr,"qim Y now %dx%d\n", qim->nx, qim->ny); */
1689 bot = (nxnew - (int)(fy*nxi))/2;
1690 zim = mri_valpad_2D( bot, nxnew-(int)(fy*nxi)-bot,
1691 0, 0, qim, pval);
1692 if (qim != zim) mri_free(qim) ;
1693 qim = zim; zim = NULL;
1694 /* fprintf(stderr,"qim Y padded %dx%d, bot=%d\n",
1695 qim->nx, qim->ny, bot); */
1696 }
1697 } else {
1698 qim = mri_resize( bim , nxnew , nynew ) ;
1699 /* fprintf(stderr,"qim now %dx%d\n", qim->nx, qim->ny); */
1700 }
1701 ADDTO_IMARR( outar , qim ) ;
1702 if( bim != imin ) mri_free(bim) ;
1703 mri_free( imin );
1704 } else {
1705 ADDTO_IMARR( outar , imin ) ;
1706 }
1707 }
1708
1709 FREE_IMARR(newar) ; /* don't need this no more */
1710 }
1711
1712 RETURN( outar );
1713 }
1714
1715 /*---------------------------------------------------------------*/
1716
1717 /*! Read a raw PPM file into 3 byte-valued MRI_IMAGEs.
1718
1719 \date 16 May 1995
1720 */
1721
mri_read_ppm3(char * fname)1722 MRI_IMARR * mri_read_ppm3( char * fname )
1723 {
1724 int ch , nch , nx,ny,maxval , length , npix,ii ;
1725 char buf[512] ;
1726 MRI_IMAGE *rim , *gim , *bim ;
1727 MRI_IMARR * outar ;
1728 FILE * imfile ;
1729 byte * rby , * gby , * bby , * rgby ;
1730
1731 ENTRY("mri_read_ppm3") ;
1732
1733 /*** open input file ***/
1734
1735 imfile = fopen( fname , "r" ) ;
1736 if( imfile == NULL ){
1737 ERROR_message("mri_read_ppm3: couldn't open file %s" , fname ) ; RETURN(NULL) ;
1738 }
1739
1740 /*** check if a raw PPM file ***/
1741
1742 ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; RETURN(NULL); }
1743 ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; RETURN(NULL); }
1744
1745 /* magic P6 found, so read numbers in header */
1746
1747 ch = getc(imfile) ;
1748
1749 NUMSCAN(nx) ; if( nx <= 0 ) { fclose(imfile) ; RETURN(NULL); }
1750 NUMSCAN(ny) ; if( ny <= 0 ) { fclose(imfile) ; RETURN(NULL); }
1751 NUMSCAN(maxval) ; if( maxval <= 0 ||
1752 maxval > 255 ) { fclose(imfile) ; RETURN(NULL); }
1753
1754 /*** create output images and workspace array ***/
1755
1756 rim = mri_new( nx , ny , MRI_byte ) ; rby = mri_data_pointer( rim ) ;
1757 gim = mri_new( nx , ny , MRI_byte ) ; gby = mri_data_pointer( gim ) ;
1758 bim = mri_new( nx , ny , MRI_byte ) ; bby = mri_data_pointer( bim ) ;
1759
1760 sprintf(buf,"%s#R",fname) ; mri_add_name( buf , rim ) ;
1761 sprintf(buf,"%s#G",fname) ; mri_add_name( buf , gim ) ;
1762 sprintf(buf,"%s#B",fname) ; mri_add_name( buf , bim ) ;
1763
1764 rgby = (byte *) malloc( sizeof(byte) * 3*nx*ny ) ;
1765 if( rgby == NULL ){
1766 fprintf(stderr,"couldn't malloc workspace in mri_read_ppm3!\n") ; EXIT(1) ;
1767 }
1768
1769 /*** read all data into workspace array ***/
1770
1771 length = fread( rgby , sizeof(byte) , 3*nx*ny , imfile ) ;
1772 fclose( imfile ) ;
1773
1774 if( length != 3*nx*ny ){
1775 free(rgby) ; mri_free(rim) ; mri_free(gim) ; mri_free(bim) ;
1776 fprintf(stderr,"couldn't read data from file %s in mri_read_ppm3\n",fname) ;
1777 RETURN(NULL);
1778 }
1779
1780 /*** put data from workspace array into output images ***/
1781
1782 npix = nx*ny ;
1783 for( ii=0 ; ii < npix ; ii++ ){
1784 rby[ii] = rgby[3*ii ] ;
1785 gby[ii] = rgby[3*ii+1] ;
1786 bby[ii] = rgby[3*ii+2] ;
1787 }
1788 free( rgby ) ;
1789
1790 /*** create output image array ***/
1791
1792 INIT_IMARR(outar) ;
1793 ADDTO_IMARR( outar , rim ) ;
1794 ADDTO_IMARR( outar , gim ) ;
1795 ADDTO_IMARR( outar , bim ) ;
1796 RETURN(outar);
1797 }
1798
1799 /*-----------------------------------------------------------------
1800 routines added 1 Oct 1995
1801 -------------------------------------------------------------------*/
1802
1803 /*! Read 1 2D image, then "nsize" it - make it a power of 2 in sizes.
1804
1805 This was developed in the days when FD/FD2/fim ruled the world, and
1806 those programs (AJ's legacy) only deal with square images that are
1807 a power of 2 in size.
1808 \date 01 Oct 1995
1809 */
1810
mri_read_nsize(char * fname)1811 MRI_IMAGE *mri_read_nsize( char * fname )
1812 {
1813 MRI_IMARR *imar ;
1814 MRI_IMAGE *imout ;
1815
1816 imar = mri_read_file( fname ) ;
1817 if( imar == NULL ) return NULL ;
1818 if( imar->num != 1 ){ DESTROY_IMARR(imar) ; return NULL ; }
1819
1820 imout = mri_nsize( IMAGE_IN_IMARR(imar,0) ) ;
1821 mri_add_name( IMAGE_IN_IMARR(imar,0)->name , imout ) ;
1822
1823 DESTROY_IMARR(imar) ;
1824 return imout ;
1825 }
1826
1827 /*! Read many 2D images from many files. */
1828
mri_read_many_nsize(int nf,char * fn[])1829 MRI_IMARR *mri_read_many_nsize( int nf , char * fn[] )
1830 {
1831 MRI_IMARR * newar , * outar ;
1832 MRI_IMAGE * im ;
1833 int ii ;
1834
1835 newar = mri_read_many_files( nf , fn ) ;
1836 if( newar == NULL ) return NULL ;
1837
1838 INIT_IMARR(outar) ;
1839 for( ii=0 ; ii < newar->num ; ii++ ){
1840 im = mri_nsize( IMAGE_IN_IMARR(newar,ii) ) ;
1841 mri_add_name( IMAGE_IN_IMARR(newar,ii)->name , im ) ;
1842 ADDTO_IMARR(outar,im) ;
1843 mri_free( IMAGE_IN_IMARR(newar,ii) ) ;
1844 }
1845 FREE_IMARR(newar) ;
1846 return outar ;
1847 }
1848
1849 /*------------------------------------------------------------------------*/
1850
1851 /*! Set up MCW_SIZE_# database for input.
1852
1853 This implements the facility for the user to define MCW_IMSIZE_1
1854 (or AFNI_IMSIZE_1) et cetera, for pre-defining a relationship between
1855 a file size in bytes and a 3D: prefix. This function is only called
1856 once to setup the table.
1857 \date 07 Nov 95
1858 */
1859
init_MCW_sizes(void)1860 void init_MCW_sizes(void)
1861 {
1862 int num , count ;
1863 char ename[32] ;
1864 char * str ;
1865
1866 if( MCW_imsize_good >= 0 ) return ;
1867
1868 MCW_imsize_good = 0 ;
1869
1870 for( num=0 ; num < MAX_MCW_IMSIZE ; num++ ){ /* look for environment string */
1871
1872 imsize[num].size = -1 ;
1873
1874 /* try to find environment variable with the num-th name */
1875
1876 sprintf( ename , "AFNI_IMSIZE_%d", num+1 ) ; str = my_getenv(ename) ;
1877 if( str == NULL ){
1878 sprintf( ename , "MCW_IMSIZE_%d" , num+1 ) ; str = my_getenv(ename) ;
1879 }
1880 if( str == NULL ){
1881 sprintf( ename, "AFNI_IMSIZE_%02d",num+1 ) ; str = my_getenv(ename) ;
1882 }
1883 if( str == NULL ){
1884 sprintf( ename, "MCW_IMSIZE_%02d" ,num+1 ) ; str = my_getenv(ename) ;
1885 }
1886 if( str == NULL ) continue ; /* no luck */
1887
1888 imsize[num].prefix = (char *) malloc( sizeof(char) * strlen(str) ) ;
1889 if( imsize[num].prefix == NULL ){
1890 fprintf(stderr,"\n*** Can't malloc in init_MCW_sizes! ***\a\n");
1891 EXIT(1) ;
1892 }
1893
1894 if( str[0] != '%' ){ /* e.g., 16096=3D:-1:0:64:64:1: */
1895
1896 imsize[num].head = -1 ;
1897 count = sscanf( str , "%d=%s" , &(imsize[num].size) , imsize[num].prefix ) ;
1898 if( count != 2 || imsize[num].size < 2 || strlen(imsize[num].prefix) < 2 ){
1899 free( imsize[num].prefix ) ;
1900 fprintf(stderr,"bad environment %s = %s\n" ,
1901 ename , str ) ;
1902 }
1903
1904 } else { /* e.g., %16096+0=3D:0:7904:64:64: */
1905
1906 count = sscanf( str+1 , "%d+%d=%s" ,
1907 &(imsize[num].size) , &(imsize[num].head) , imsize[num].prefix ) ;
1908
1909 if( count != 3 || imsize[num].size < 2 ||
1910 imsize[num].head < 0 || strlen(imsize[num].prefix) < 2 ){
1911
1912 free( imsize[num].prefix ) ;
1913 fprintf(stderr,"bad environment %s = %s\n" ,
1914 ename , str ) ;
1915 }
1916 }
1917
1918 MCW_imsize_good ++ ;
1919 }
1920
1921 return ;
1922 }
1923
1924 /*------------------------------------------------------------------------------*/
1925 /*! My version of strdup(), which won't fail if the input is NULL. */
1926
my_strdup(char * str)1927 char * my_strdup( char * str )
1928 {
1929 char * new_str ;
1930 if( str == NULL ) return NULL ;
1931 new_str = (char *) malloc( sizeof(char) * (strlen(str)+1) ) ;
1932 if( new_str != NULL ) strcpy( new_str , str ) ;
1933 return new_str ;
1934 }
1935
1936 /*------------------------------------------------------------------------------*/
1937
1938 /*! Check if a filesize fits an MCW_IMSIZE setup.
1939
1940 \param fname = Filename
1941 \return A new "filename" with 3D header attached if it fits.
1942 If not, return a copy of the filename. In any case the
1943 returned string should be free()-d when it is no longer needed.
1944 */
1945
imsized_fname(char * fname)1946 char * imsized_fname( char * fname )
1947 {
1948 int num , lll ;
1949 long long len ; /* 22 Mar 2007 */
1950 char * new_name ;
1951
1952 init_MCW_sizes() ;
1953 if( MCW_imsize_good == 0 ){
1954 new_name = my_strdup(fname) ; /* nothing to fit */
1955 return new_name ; /* --> return copy of old name */
1956 }
1957
1958 len = THD_filesize( fname ) ;
1959 if( len <= 0 ){
1960 new_name = my_strdup(fname) ; /* not an existing filename */
1961 return new_name ; /* --> return copy of old name */
1962 }
1963
1964 for( num=0 ; num < MAX_MCW_IMSIZE ; num++ ){ /* check each possibility */
1965
1966 if( imsize[num].size <= 0 ) continue ; /* skip to next one */
1967
1968 if( imsize[num].head < 0 && len == imsize[num].size ){ /* fixed size fit */
1969
1970 lll = strlen(fname) + strlen(imsize[num].prefix) + 4 ;
1971 new_name = (char *) malloc( sizeof(char) * lll ) ;
1972 if( new_name == NULL ){
1973 fprintf(stderr,"\n*** Can't malloc in imsized_fname! ***\a\n");
1974 EXIT(1) ;
1975 }
1976 sprintf( new_name , "%s%s" , imsize[num].prefix , fname ) ;
1977 return new_name ;
1978
1979 } else if( (len-imsize[num].head) % imsize[num].size == 0 ){
1980 int count = (len-imsize[num].head) / imsize[num].size ;
1981
1982 if( count < 1 ) continue ; /* skip to next one */
1983
1984 lll = strlen(fname) + strlen(imsize[num].prefix) + 32 ;
1985 new_name = (char *) malloc( sizeof(char) * lll ) ;
1986 if( new_name == NULL ){
1987 fprintf(stderr,"\n*** Can't malloc in imsized_fname! ***\a\n");
1988 EXIT(1) ;
1989 }
1990 sprintf( new_name , "%s%d:%s" , imsize[num].prefix , count , fname ) ;
1991 return new_name ;
1992 }
1993
1994 }
1995
1996 new_name = my_strdup(fname) ; /* no fit --> return copy of old name */
1997 return new_name ;
1998 }
1999
2000 #if 0 /* removed on 22 Mar 2007 */
2001 /*------------------------------------------------------------------------*/
2002 /*! Return the size of a file in bytes.
2003
2004 \param pathname = input filename
2005 \return File length if file exists; -1 if it doesn't.
2006 \see THD_filesize() in thd_filestuff.c.
2007 */
2008
2009 long mri_filesize( char * pathname )
2010 {
2011 static struct stat buf ;
2012 int ii ;
2013
2014 if( pathname == NULL ) return -1 ;
2015 ii = stat( pathname , &buf ) ; if( ii != 0 ) return -1 ;
2016 return buf.st_size ;
2017 }
2018 #endif
2019
2020 /*---------------------------------------------------------------*/
2021
2022 /*! Read the header from PPM file and return its info.
2023
2024 \param fname = file name
2025 \return *nx and *ny are set to the image dimensions;
2026 if they are set to 0, something bad happened
2027 (e.g., the file isn't a PPM file, or doesn't exist).
2028 \date 17 Sep 2001
2029 */
2030
mri_read_ppm_header(char * fname,int * nx,int * ny)2031 void mri_read_ppm_header( char *fname , int *nx, int *ny )
2032 {
2033 FILE *imfile ;
2034 int ch , nch , nxx,nyy ;
2035 char buf[256] ;
2036
2037 ENTRY("mri_read_ppm_header") ;
2038
2039 if( fname == NULL || nx == NULL || ny == NULL ) EXRETURN ;
2040
2041 *nx = *ny = 0 ; /* default returns */
2042
2043 /*** open input file ***/
2044
2045 imfile = fopen( fname , "r" ) ;
2046 if( imfile == NULL ){
2047 ERROR_message("mri_read_ppm_header: couldn't open file %s" , fname ) ; EXRETURN ;
2048 }
2049
2050 /*** check if a raw PPM file ***/
2051
2052 ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; EXRETURN ; }
2053 ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; EXRETURN ; }
2054
2055 /* magic P6 found, so read numbers in header */
2056
2057 ch = getc(imfile) ;
2058
2059 NUMSCAN(nxx) ; if( nxx <= 0 ){ fclose(imfile) ; EXRETURN ; }
2060 NUMSCAN(nyy) ; if( nyy <= 0 ){ fclose(imfile) ; EXRETURN ; }
2061
2062 /* return dimensions */
2063
2064 fclose(imfile) ; *nx = nxx ; *ny = nyy ; EXRETURN ;
2065 }
2066
2067 /*---------------------------------------------------------------*/
2068
2069 /*! Reads a raw PPM file into 1 2D MRI_rgb-valued image.
2070
2071 \param fname = Image filename
2072 \return An MRI_IMAGE if things worked OK; NULL if not
2073 \date 13 May 1996
2074 */
2075
mri_read_ppm(char * fname)2076 MRI_IMAGE * mri_read_ppm( char * fname )
2077 {
2078 int ch , nch , nx,ny,maxval , length ;
2079 MRI_IMAGE * rgbim ;
2080 FILE * imfile ;
2081 byte * rgby ;
2082 char buf[256] ;
2083
2084 ENTRY("mri_read_ppm") ;
2085
2086 /*** open input file ***/
2087
2088 imfile = fopen( fname , "r" ) ;
2089 if( imfile == NULL ){
2090 ERROR_message("mri_read_ppm: couldn't open file %s" , fname ) ; RETURN(NULL) ;
2091 }
2092
2093 /*** check if a raw PPM file ***/
2094
2095 ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; RETURN(NULL); }
2096 ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; RETURN(NULL); }
2097
2098 /* magic P6 found, so read numbers in header */
2099
2100 ch = getc(imfile) ;
2101
2102 NUMSCAN(nx) ; if( nx <= 0 ) { fclose(imfile); RETURN(NULL); }
2103 NUMSCAN(ny) ; if( ny <= 0 ) { fclose(imfile); RETURN(NULL); }
2104 NUMSCAN(maxval); if( maxval <= 0 ||
2105 maxval > 255 ){ fclose(imfile); RETURN(NULL); }
2106
2107 /*** create output image ***/
2108
2109 rgbim = mri_new( nx , ny , MRI_rgb ) ; mri_add_name( fname , rgbim ) ;
2110 rgby = MRI_RGB_PTR(rgbim) ;
2111
2112 /*** read all data into image array */
2113
2114 length = fread( rgby , sizeof(byte) , 3*nx*ny , imfile ) ;
2115 fclose( imfile ) ;
2116
2117 if( length != 3*nx*ny ){ mri_free(rgbim) ; RETURN(NULL) ; }
2118
2119 /* 17 Sep 2001: scale to maxval=255, if needed */
2120
2121 if( maxval < 255 ){
2122 int ii ; float fac = 255.4/maxval ;
2123 for( ii=0 ; ii < 3*nx*ny ; ii++ ) rgby[ii] = (byte)( rgby[ii]*fac ) ;
2124 }
2125
2126 RETURN(rgbim) ;
2127 }
2128 /*---------------------------------------------------------------*/
2129
2130 /*! Length of line buffer for mri_read_ascii() */
2131 /* rcr - improve this */
2132 #define LBUF 5048576 /* 08 Jul 2004: increased to 512K from 64K
2133 27 Dec 2012: increased to 1024K from 512K */
2134
2135 /*! Free a buffer and set it to NULL */
2136 #define FRB(b) do{ if( (b)!=NULL ){free((b)); (b)=NULL;} }while(0)
2137
2138 #undef USE_LASTBUF
2139
2140 static int save_comments = 0 ; /* 03 Aug 2016 */
2141 static char *comment_buffer = NULL ;
2142
2143 /*---------------------------------------------------------------*/
2144 /*! [20 Jun 2002] Like fgets, but also
2145 - skips blank or comment lines
2146 - skips leading and trailing whitespace
2147 - catenates lines that end in '\' (replacing '\' with ' ')
2148 - returns duplicate of last line if first 2
2149 nonblank input characters are "" [20 Jul 2004]
2150 -----------------------------------------------------------------*/
2151
my_fgets(char * buf,int size,FILE * fts)2152 static char * my_fgets( char *buf , int size , FILE *fts )
2153 {
2154 char *ptr ;
2155 int nbuf , ll,ii , cflag ;
2156 static char *qbuf=NULL ;
2157
2158 #ifdef USE_LASTBUF
2159 static char *lastbuf = NULL ; /* 20 Jul 2004 */
2160 static int nlastbuf = 0 ;
2161
2162 if( buf == NULL && lastbuf != NULL ){ /* 20 Jul 2004 */
2163 free((void *)lastbuf); lastbuf = NULL; nlastbuf = 0 ;
2164 }
2165 #endif
2166
2167 if( buf == NULL && qbuf != NULL ){ free((void *)qbuf); qbuf = NULL; }
2168
2169 if( buf == NULL || size < 1 || fts == NULL ){
2170 FRB(comment_buffer) ; return NULL ;
2171 }
2172
2173 if( qbuf == NULL ) qbuf = AFMALL(char, LBUF) ; /* 1st time in */
2174
2175 nbuf = 0 ; /* num bytes stored in buf so far */
2176 cflag = 0 ; /* flag if we're catenating lines */
2177
2178 while(1){ /* loop and read lines, creating a logical line */
2179
2180 ptr = afni_fgets( qbuf , LBUF , fts ) ; /* read next whole line */
2181
2182 if( ptr == NULL ) break ; /* must be end-of-file */
2183
2184 /* skip leading whitespace */
2185
2186 for( ; *ptr != '\0' && isspace(*ptr) ; ptr++ ) ; /* nada */
2187
2188 /* skip entirely blank lines, unless we are catenating */
2189
2190 if( *ptr == '\0' ){ if(cflag) break; else continue; }
2191
2192 #ifdef USE_LASTBUF
2193 /* if a duplicate is requested, return it now [20 Jul 2004] */
2194
2195 if( *ptr == '"' && *(ptr+1) == '"' && nlastbuf > 0 && nbuf == 0 ){
2196 ll = strlen(lastbuf) ; if( ll >= size ) ll = size-1 ;
2197 memcpy(buf,lastbuf,ll-1) ; buf[ll] = '\0' ;
2198 return buf ;
2199 }
2200 #endif
2201
2202 /* strip trailing whitespace */
2203
2204 ll = strlen(ptr) ; /* will be > 0 */
2205 for( ii=ll-1 ; isspace(ptr[ii]) && ii > 0 ; ii-- ) /* blank => NUL */
2206 ptr[ii] = '\0' ;
2207
2208 ll = strlen(ptr) ; /* number of chars left */
2209 if( ll == 0 ) continue ; /* should not happen */
2210
2211 /* save then skip comment lines (even if we are catenating) */
2212
2213 if( *ptr == '#' || (*ptr == '/' && *(ptr+1) == '/') ){
2214 if( save_comments ){
2215 comment_buffer = THD_zzprintf( comment_buffer , "%s\n" , ptr ) ;
2216 }
2217 continue ;
2218 }
2219
2220 cflag = (ptr[ll-1] == '\\') ; /* catenate next line? */
2221 if( cflag ) ptr[ll-1] = ' ' ; /* replace '\' with ' ' */
2222
2223 /* now copy what's left (ll+1 bytes) at tail of output buffer */
2224
2225 if( nbuf+ll+1 > size ){ /* too much for output buffer? */
2226 ll = size - (nbuf+1) ;
2227 if( ll <= 0 ) break ; /* should not happen */
2228 }
2229
2230 memcpy(buf+nbuf,ptr,ll+1) ; nbuf += ll ;
2231 if( !cflag ) break ;
2232
2233 } /* loop to get next line if catenation is turned on */
2234
2235 #ifdef LASTBUF
2236 /* make a copy of result in lastbuf [20 Jul 2004] */
2237
2238 ll = strlen(buf) ;
2239 if( ll+1 > nlastbuf ){
2240 nlastbuf = ll+2 ; lastbuf = (char *)realloc((void *)lastbuf,nlastbuf) ;
2241 }
2242 memcpy(lastbuf,buf,ll+1) ;
2243 #endif
2244
2245 /* and we is done */
2246
2247 if( nbuf > 0 ) return buf ; /* return what we read already */
2248 return NULL ; /* signal of failure get data */
2249 }
2250
2251 /*--------------------------------------------------------------*/
2252 static float lbfill = 0.0f ; /* 10 Aug 2004 */
2253 static int oktext = 0; /* ZSS: Oct 16 2009 */
2254 static int linebufdied = 0; /* ZSS: Oct 19 2009 */
2255 static int doublelinebufdied = 0;/* ZSS: Oct 19 2009 */
2256 /*--------------------------------------------------------------*/
2257
2258 /* Consider an 'i' to be a complex number flag if it is right after a digit
2259 ZSS Jan 2011*/
2260 #define ISCOMPLEXi(b,t) ((b[t] == 'i' && t != 0 && b[t-1]>= '0' && b[t-1]<='9'))
2261
2262 /* return a 1 if c is a not a valid first non-white char of a
2263 non-comment 1D line */
iznogood_1D(char * cv,int t)2264 byte iznogood_1D (char *cv, int t)
2265 {
2266 char c=cv[t];
2267 if ( (c < '0' || c > '9') &&
2268 c != '+' && c != '-' && c != '.' && c != 'e' &&
2269 !ISCOMPLEXi(cv,t) && c != ',' && /* allow for complex input */
2270 c != '@' && c != '*' /* allow for special 1D trickery */
2271 ) return 1;
2272 else return 0;
2273
2274 }
2275
2276 /*! Decode a line buffer into an array of floats. */
decode_linebuf(char * buf)2277 static floatvec * decode_linebuf( char *buf ) /* 20 Jul 2004 */
2278 {
2279 floatvec *fv=NULL ;
2280 int blen, bpos, ncol, ii, jj, temppos, count;
2281 int alloc_chunk, alloc_unit = 10000, incr;
2282 int n_alloced = 0, slowmo = 0 ; /* ZSS speedups */
2283 char sep, vbuf[64] , *cpt, *ope=NULL;
2284 float val ;
2285
2286 if( buf == NULL || *buf == '\0' ) return fv ;
2287
2288 blen = strlen(buf) ;
2289 ncol = 0 ;
2290 linebufdied = 0;
2291
2292 /* convert commas (or 'i' for complex numbers ZSS Oct 06) to blanks */
2293 /* note 'e' is commonly found in numeric files as in scientific notation*/
2294 for( ii=0 ; ii < blen ; ii++ ) {
2295 temppos = ii; incr = 0;
2296 if(isalpha(buf[ii])){
2297 /* skip past alphabetics in a row*/
2298 jj = ii;
2299 for( ; jj < blen && isalpha(buf[jj]) ; jj++ ) ;
2300 incr = jj - ii - 1; /* only move if more than 1 char long */
2301 if(incr) ii = jj;
2302 }
2303
2304 /* convert some alphabetic characters to space (:,i)
2305 if they are not followed by other alphabetics */
2306 if( incr<=0 &&
2307 ( buf[temppos] == ',' || ISCOMPLEXi(buf, temppos)
2308 || buf[temppos] == ':' ) )
2309 buf[temppos] = ' ' ;
2310 /* turn on "slow mo" reading if non-numeric */
2311 if( !slowmo &&
2312 (buf[temppos] == '*' || buf[temppos] == '@' ||
2313 isalpha(buf[temppos])) ) slowmo = 1;
2314 }
2315
2316 fv = (floatvec *)malloc(sizeof(floatvec)) ;
2317 fv->nar = 0 ;
2318 fv->ar = (float *)NULL ;
2319
2320 for( bpos=0 ; bpos < blen ; ){
2321 /* skip to next nonblank character */
2322
2323 for( ; bpos < blen && isspace(buf[bpos]) ; bpos++ ) ; /* nada */
2324 if( bpos == blen ) break ; /* end of line */
2325
2326
2327 val = 0.0 ; count = 1 ;
2328 if (slowmo) { /* trickery */
2329 sscanf( buf+bpos , "%63s" , vbuf ) ;
2330 if (!oktext && iznogood_1D(buf, bpos)) {/* Morality Police Oct 16 09 */
2331 if (vbuf[0] != '#') { /* not a comment, die */
2332 linebufdied = 1;
2333 fv->nar = 0; /* this will cause a clean up on the way out */
2334 /* By setting nar to 0, reading of 1D file will terminate*/
2335 /* at first row where text is encountered whenever */
2336 /* AFNI_1D_ZERO_TEXT is not YES . This is more consistent */
2337 /* than the alternative. For example: */
2338 /* if nar is not set to 0, then a file that has */
2339 /* 1 a 3
2340 4 b 6 comes out as a 1 column vector,
2341 but
2342 1 2 3
2343 4 b 6 comes out as a 3x3 matrix */
2344 }
2345 break;
2346 }
2347
2348 if( vbuf[0] == '*' || isalpha(vbuf[0]) ){ /* 10 Aug 2004 */
2349 val = lbfill ;
2350 } else if( (cpt=strchr(vbuf,'@')) != NULL ){
2351 sscanf( vbuf , "%d%c%f" , &count , &sep , &val ) ;
2352 if( count < 1 ) count = 1 ;
2353 if( *(cpt+1) == '*' ) val = lbfill ; /* 10 Aug 2004 */
2354 } else {
2355 sscanf( vbuf , "%f" , &val ) ;
2356 }
2357 incr = strlen(vbuf) ;
2358 } else { /* no muss no fuss, take it straight */
2359 /* sscanf( vbuf , "%f" , &val ) ; slow, slow, tan go close*/
2360 val = strtod(buf+bpos, &ope);
2361 incr = ope - (buf+bpos);
2362 }
2363 if( incr <= 0 ) break ; /* 16 Oct 2007 */
2364 if (fv->nar+count > n_alloced) {
2365 /* fprintf(stderr,"reallocing past %d with count %d...\n", n_alloced, count); */
2366 if (count > alloc_unit) alloc_chunk = count;
2367 else alloc_chunk = alloc_unit;
2368 fv->ar = (float *)realloc( (void *)fv->ar , sizeof(float)*(n_alloced+alloc_chunk) );
2369 n_alloced = n_alloced + alloc_chunk;
2370 }
2371 for( ii=0 ; ii < count ; ii++ ) fv->ar[ii+fv->nar] = val ;
2372 fv->nar += count ;
2373 bpos += incr ;
2374 }
2375
2376 if( fv->nar == 0 ){ KILL_floatvec(fv); fv = NULL; }
2377 else { if (fv->nar < n_alloced) fv->ar = (float *)realloc( (void *)fv->ar , sizeof(float)*(fv->nar) ); }
2378 return fv ;
2379 }
2380
decode_double_linebuf(char * buf)2381 static doublevec * decode_double_linebuf( char *buf ) /* 20 Jul 2004 */
2382 {
2383 doublevec *dv=NULL ;
2384 int blen, bpos, ncol, ii, jj, temppos, count ;
2385 int alloc_chunk, alloc_unit = 10000, incr;
2386 int n_alloced = 0, slowmo = 0 ; /* ZSS speedups */
2387 char sep, vbuf[64] , *cpt , *ope=NULL;
2388 double val ;
2389
2390 if( buf == NULL || *buf == '\0' ) return dv ;
2391
2392 blen = strlen(buf) ;
2393 ncol = 0 ;
2394 doublelinebufdied = 0;
2395
2396 /* convert commas (or 'i' for complex numbers ZSS Oct 06) to blanks */
2397 /* note 'e' is commonly found in numeric files as in scientific notation*/
2398 for( ii=0 ; ii < blen ; ii++ ) {
2399 temppos = ii; incr = 0;
2400 if(isalpha(buf[ii])){
2401 /* skip past alphabetics in a row*/
2402 jj = ii;
2403 for( ; jj < blen && isalpha(buf[jj]) ; jj++ ) ;
2404 incr = jj - ii - 1; /* only move if more than 1 char long */
2405 if(incr) ii = jj;
2406 }
2407
2408 /* convert some alphabetic characters to space (:,i)
2409 if they are not followed by other alphabetics */
2410 if( incr<=0 &&
2411 ( buf[temppos] == ',' || ISCOMPLEXi(buf, temppos)
2412 || buf[temppos] == ':' ) )
2413 buf[temppos] = ' ' ;
2414 /* turn on "slow mo" reading if non-numeric */
2415 if( !slowmo &&
2416 (buf[temppos] == '*' || buf[temppos] == '@' ||
2417 isalpha(buf[temppos])) ) slowmo = 1;
2418 }
2419
2420 dv = (doublevec *)malloc(sizeof(doublevec)) ;
2421 dv->nar = 0 ;
2422 dv->ar = (double *)NULL ;
2423
2424 for( bpos=0 ; bpos < blen ; ){
2425 /* skip to next nonblank character */
2426
2427 for( ; bpos < blen && isspace(buf[bpos]) ; bpos++ ) ; /* nada */
2428 if( bpos == blen ) break ; /* end of line */
2429
2430
2431 val = 0.0 ; count = 1 ;
2432 if (slowmo) { /* trickery */
2433 sscanf( buf+bpos , "%63s" , vbuf ) ;
2434 if (!oktext && iznogood_1D(buf, bpos)) {/* Morality Police Oct 16 09 */
2435 if (vbuf[0] != '#') { /* not a comment, die */
2436 doublelinebufdied = 1;
2437 dv->nar = 0; /* for comment, see same section in decode_linebuf */
2438 }
2439 break;
2440 }
2441
2442 if( vbuf[0] == '*' || isalpha(vbuf[0]) ){ /* 10 Aug 2004 */
2443 val = (double)lbfill ;
2444 } else if( (cpt=strchr(vbuf,'@')) != NULL ){
2445 sscanf( vbuf , "%d%c%lf" , &count , &sep , &val ) ;
2446 if( count < 1 ) count = 1 ;
2447 if( *(cpt+1) == '*' ) val = (double)lbfill ; /* 10 Aug 2004 */
2448 } else {
2449 sscanf( vbuf , "%lf" , &val ) ;
2450 }
2451 incr = strlen(vbuf) ;
2452 } else { /* no muss no fuss, take it straight */
2453 /* sscanf( vbuf , "%f" , &val ) ; slow, slow, tan go close*/
2454 val = strtod(buf+bpos, &ope);
2455 incr = ope - (buf+bpos);
2456 }
2457 if (dv->nar+count > n_alloced) {
2458 /* fprintf(stderr,"reallocing past %d with count %d...\n", n_alloced, count); */
2459 if (count > alloc_unit) alloc_chunk = count;
2460 else alloc_chunk = alloc_unit;
2461 dv->ar = (double *)realloc( (void *)dv->ar , sizeof(double)*(n_alloced+alloc_chunk) );
2462 n_alloced = n_alloced + alloc_chunk;
2463 }
2464
2465 for( ii=0 ; ii < count ; ii++ ) dv->ar[ii+dv->nar] = val ;
2466 dv->nar += count ;
2467 bpos += incr ;
2468 }
2469
2470 if( dv->nar == 0 ){ KILL_doublevec(dv); dv = NULL; }
2471 else { if (dv->nar < n_alloced) dv->ar = (double *)realloc( (void *)dv->ar , sizeof(double)*(dv->nar) ); }
2472 return dv ;
2473 }
2474
2475 /*---------------------------------------------------------------*/
2476
2477 /*! Increment for time series array size for mri_read_ascii() */
2478 #define INC_TSARSIZE 128
2479
2480 /*! Read an array of ASCII numbers into a 1D or 2D image.
2481
2482 \param fname = input filename
2483 \return Pointer to MRI_IMAGE (in MRI_float) format if things
2484 are cool; NULL if not.
2485 \date Jun 1996
2486
2487 Example input:
2488 - Line 1: 3 4 6
2489 - Line 2: 2 2 2
2490 - Line 3: 7 2 1
2491 - Line 4: 9 9 6
2492 This produces an image with nx=3 and ny=4. The first row
2493 is read to determine nx; all subsequent rows must have nx
2494 values. A line whose very first character is a '#' will
2495 be skipped as a comment. A line with no characters (just
2496 the '\n') will also be skipped.
2497
2498 20 Jun 2002: modified to use my_fgets() instead of fgets().
2499 */
2500
mri_read_ascii(char * fname)2501 MRI_IMAGE * mri_read_ascii( char *fname )
2502 {
2503 MRI_IMAGE *outim ;
2504 int ii,jj,val , used_tsar , alloc_tsar ;
2505 float *tsar ;
2506 float ftemp ;
2507 FILE *fts ;
2508 char *ptr ;
2509 int ncol , bpos , blen , nrow ;
2510 static char *buf=NULL ; /* 20 Jun 2002: make a ptr */
2511
2512 floatvec *fvec ; /* 20 Jul 2004 */
2513 int incts ;
2514 char *cbuf=NULL ; /* 03 Aug 2016: comments? */
2515
2516 ENTRY("mri_read_ascii") ;
2517
2518 if (AFNI_yesenv("AFNI_1D_ZERO_TEXT")) oktext = 1; /* ZSS Oct 16 09 */
2519 else oktext = 0;
2520
2521 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
2522
2523 STATUS(fname) ; /* 16 Oct 2007 */
2524
2525 if( strncmp(fname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
2526 MRI_IMAGE *qim = mri_1D_fromstring( fname+3 ) ;
2527 if( qim != NULL ){
2528 outim = mri_transpose(qim); mri_free(qim); RETURN(outim);
2529 }
2530 }
2531
2532 fts = fopen( fname , "r" );
2533 if( fts == NULL ){ NI_sleep(33); fts = fopen( fname , "r" ); }
2534 if( fts == NULL ){
2535 ERROR_message("mri_read_ascii: couldn't open file %s" , fname ) ; RETURN(NULL) ;
2536 }
2537
2538 if( buf == NULL ) buf = AFMALL(char, LBUF) ; /* create buffer */
2539
2540 /** step 1: read in the first line and see how many numbers are in it
2541 (skipping lines that are comments or entirely blank) */
2542
2543 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2544
2545 save_comments = 1 ; /* 03 Aug 2016 */
2546 ptr = my_fgets( buf , LBUF , fts ) ;
2547 if( ptr==NULL || *ptr=='\0' ){ /* bad read? */
2548 FRB(comment_buffer); FRB(buf); fclose(fts);
2549 ERROR_message("mri_read_ascii: can't read any valid data from file %s",fname) ;
2550 RETURN(NULL) ;
2551 }
2552 save_comments = 0;
2553 if( comment_buffer ) cbuf = strdup(comment_buffer); /* 6 Aug 2016 [rickr] */
2554 else cbuf = NULL;
2555 FRB(comment_buffer);
2556
2557 lbfill = 0.0f ; /* 10 Aug 2004 */
2558
2559 fvec = decode_linebuf( buf ) ; /* 20 Jul 2004 */
2560 if( fvec == NULL || fvec->nar == 0 ){
2561 if (strlen(buf) > 64) {
2562 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2563 }
2564 if (linebufdied) {/* death? ZSS: Oct 19 2009*/
2565 static int nfail=0 ;
2566 if( ++nfail < 9 ){
2567 fprintf(stderr,
2568 "\n** Error: Failed parsing data row 0 of 1D file '%.44s'\n"
2569 " Check for illegal non-numeric characters in:\n"
2570 " '%s'\n", fname,buf);
2571 if (nfail == 1)
2572 fprintf(stderr,
2573 " You can have text entries set to 0 with -ok_1D_text or by \n"
2574 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2575
2576 }
2577 }
2578 if( fvec != NULL ) KILL_floatvec(fvec) ;
2579 FRB(buf); fclose(fts); RETURN(NULL);
2580 }
2581 ncol = fvec->nar ; KILL_floatvec(fvec) ;
2582
2583 /** At this point, ncol is the number of floats to be read from each line **/
2584
2585 rewind( fts ) ; /* will start over */
2586
2587 incts = MAX(INC_TSARSIZE,ncol) ;
2588 used_tsar = 0 ;
2589 alloc_tsar = incts ;
2590 tsar = (float *) malloc( sizeof(float) * alloc_tsar ) ;
2591 if( tsar == NULL ){
2592 fprintf(stderr,"\n*** malloc error in mri_read_ascii ***\n"); EXIT(1);
2593 }
2594
2595 /** read lines, convert to floats, store **/
2596
2597 nrow = 0 ;
2598 while( 1 ){
2599 ptr = my_fgets( buf , LBUF , fts ) ; /* read, skipping comments*/
2600 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
2601
2602 fvec = decode_linebuf( buf ) ;
2603 if( fvec == NULL ) break ;
2604 if( fvec->nar == 0 ){ KILL_floatvec(fvec); break; }
2605
2606 if( used_tsar + ncol >= alloc_tsar ){
2607 alloc_tsar += incts ;
2608 tsar = (float *)realloc( (void *)tsar,sizeof(float)*alloc_tsar );
2609 if( tsar == NULL ){
2610 fprintf(stderr,"\n*** realloc error in mri_read_ascii ***\n"); EXIT(1);
2611 }
2612 }
2613 for( ii=0 ; ii < fvec->nar && ii < ncol ; ii++ )
2614 tsar[used_tsar+ii] = fvec->ar[ii] ;
2615 for( ; ii < ncol ; ii++ )
2616 tsar[used_tsar+ii] = 0.0 ;
2617 used_tsar += ncol ;
2618 KILL_floatvec(fvec) ;
2619
2620 nrow++ ; /* got one more complete row! */
2621 }
2622 fclose( fts ) ; /* finished with this file! */
2623 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2624
2625 if (linebufdied) {/* death? ZSS: Oct 19 2009 */
2626 static int nfail=0 ;
2627 if (strlen(buf) > 64) {
2628 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2629 }
2630 if( ++nfail < 9 ) {
2631 fprintf(stderr,
2632 "\n** Error: Failed parsing data row %d of 1D file '%.44s'\n"
2633 " Check for illegal non-numeric characters in:\n"
2634 " '%s'\n", nrow, fname , buf);
2635 if (nfail == 1)
2636 fprintf(stderr,
2637 " You can have text entries set to 0 with -ok_1D_text or by \n"
2638 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2639 }
2640 if (tsar) free(tsar); tsar = NULL;
2641 FRB(buf);
2642 RETURN(NULL);
2643 }
2644
2645 /* from <= 1 to < 1 (allow 1x1 image) 25 Jan 2006 [rickr] */
2646 if( used_tsar < 1 ){ FRB(buf); free(tsar); RETURN(NULL); }
2647
2648 tsar = (float *) realloc( tsar , sizeof(float) * used_tsar ) ;
2649 if( tsar == NULL ){
2650 fprintf(stderr,"\n*** final realloc error in mri_read_ascii ***\n"); EXIT(1);
2651 }
2652
2653 outim = mri_new_vol_empty( ncol , nrow , 1 , MRI_float ) ;
2654 mri_fix_data_pointer( tsar , outim ) ;
2655 mri_add_name( fname , outim ) ;
2656
2657 if( cbuf != NULL ){ /* 03 Aug 2016 */
2658 outim->comments = cbuf ; cbuf = NULL ;
2659 }
2660
2661 FRB(buf) ; RETURN(outim) ;
2662 }
2663
2664 /*---------------------------------------------------------------*/
2665
mri_read_double_ascii(char * fname)2666 MRI_IMAGE * mri_read_double_ascii( char * fname )
2667 {
2668 MRI_IMAGE * outim ;
2669 int ii,jj,val , used_tsar , alloc_tsar ;
2670 double * dtsar ;
2671 double dtemp ;
2672 FILE * fts ;
2673 char * ptr ;
2674 int ncol , bpos , blen , nrow ;
2675 static char *buf=NULL ; /* 20 Jun 2002: make a ptr */
2676
2677 doublevec *dvec ; /* 20 Jul 2004 */
2678 int incts ;
2679
2680 ENTRY("mri_read_double_ascii") ;
2681
2682 if (AFNI_yesenv("AFNI_1D_ZERO_TEXT")) oktext = 1; /* ZSS Oct 16 09 */
2683 else oktext = 0;
2684
2685 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
2686
2687 if( strncmp(fname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
2688 /*
2689 MRI_IMAGE *qim = mri_1D_double_fromstring( fname+3 ) ;
2690 if( qim != NULL ){
2691 outim = mri_transpose(qim); mri_free(qim); RETURN(outim);
2692 }*/
2693 fprintf(stderr,"Somebody was too lazy to allow this option here.\n"); RETURN(NULL);
2694 }
2695
2696 fts = fopen( fname , "r" );
2697 if( fts == NULL ){
2698 ERROR_message("mri_read_double_ascii: couldn't open file %s" , fname ) ; RETURN(NULL) ;
2699 }
2700
2701 if( buf == NULL ) buf = AFMALL(char, LBUF) ; /* create buffer */
2702
2703 /** step 1: read in the first line and see how many numbers are in it
2704 (skipping lines that are comments or entirely blank) */
2705
2706 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2707 ptr = my_fgets( buf , LBUF , fts ) ;
2708 if( ptr==NULL || *ptr=='\0' ){ FRB(buf); fclose(fts); RETURN(NULL); } /* bad read? */
2709
2710 lbfill = 0.0f ; /* 10 Aug 2004 */
2711
2712 dvec = decode_double_linebuf( buf ) ; /* 20 Jul 2004 */
2713 if( dvec == NULL || dvec->nar == 0 ){
2714 if (strlen(buf) > 64) {
2715 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2716 }
2717 if (doublelinebufdied) {/* death? ZSS: Oct 19 2009 */
2718 static int nfail=0 ;
2719 if( ++nfail < 9 ) {
2720 fprintf(stderr,
2721 "\n** Error: Failed parsing data row 0 of 1D file '%.44s'\n"
2722 " Check for illegal non-numeric characters in:\n"
2723 " '%s'\n", fname,buf);
2724 if (nfail == 1)
2725 fprintf(stderr,
2726 " You can have text entries set to 0 with -ok_1D_text or by \n"
2727 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2728 }
2729 }
2730 if( dvec != NULL ) KILL_doublevec(dvec) ;
2731 FRB(buf); fclose(fts); RETURN(NULL);
2732 }
2733 ncol = dvec->nar ; KILL_doublevec(dvec) ;
2734
2735 /** At this point, ncol is the number of floats to be read from each line **/
2736
2737 rewind( fts ) ; /* will start over */
2738
2739 incts = MAX(INC_TSARSIZE,ncol) ;
2740 used_tsar = 0 ;
2741 alloc_tsar = incts ;
2742 dtsar = (double *) malloc( sizeof(double) * alloc_tsar ) ;
2743 if( dtsar == NULL ){
2744 fprintf(stderr,"\n*** malloc error in mri_read_double_ascii ***\n"); EXIT(1);
2745 }
2746
2747 /** read lines, convert to floats, store **/
2748
2749 nrow = 0 ;
2750 while( 1 ){
2751 ptr = my_fgets( buf , LBUF , fts ) ; /* read */
2752 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
2753
2754 dvec = decode_double_linebuf( buf ) ;
2755 if( dvec == NULL ) break ;
2756 if( dvec->nar == 0 ){ KILL_doublevec(dvec); break; }
2757
2758 if( used_tsar + ncol >= alloc_tsar ){
2759 alloc_tsar += incts ;
2760 dtsar = (double *)realloc( (void *)dtsar,sizeof(double)*alloc_tsar );
2761 if( dtsar == NULL ){
2762 fprintf(stderr,"\n*** realloc error in mri_read_double_ascii ***\n"); EXIT(1);
2763 }
2764 }
2765 for( ii=0 ; ii < dvec->nar && ii < ncol ; ii++ )
2766 dtsar[used_tsar+ii] = dvec->ar[ii] ;
2767 for( ; ii < ncol ; ii++ )
2768 dtsar[used_tsar+ii] = 0.0 ;
2769 used_tsar += ncol ;
2770 KILL_doublevec(dvec) ;
2771
2772 nrow++ ; /* got one more complete row! */
2773 }
2774 fclose( fts ) ; /* finished with this file! */
2775 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2776
2777 if (doublelinebufdied) {/* death? ZSS: Oct 19 2009 */
2778 static int nfail=0 ;
2779 if (strlen(buf) > 64) {
2780 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2781 }
2782 if( ++nfail < 9 ) {
2783 fprintf(stderr,
2784 "\n** Error: Failed parsing data row %d of 1D file '%.44s'\n"
2785 " Check for illegal non-numeric characters in:\n"
2786 " '%s'\n", nrow, fname , buf);
2787 if (nfail == 1)
2788 fprintf(stderr,
2789 " You can have text entries set to 0 with -ok_1D_text or by \n"
2790 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2791 }
2792 if (dtsar) free(dtsar); dtsar = NULL;
2793 FRB(buf);
2794 RETURN(NULL);
2795 }
2796
2797 /* from <= 1 to < 1 (allow 1x1 image) 25 Jan 2006 [rickr] */
2798 if( used_tsar < 1 ){ FRB(buf); free(dtsar); RETURN(NULL); }
2799
2800 dtsar = (double *) realloc( dtsar , sizeof(double) * used_tsar ) ;
2801 if( dtsar == NULL ){
2802 fprintf(stderr,"\n*** final realloc error in mri_read_double_ascii ***\n"); EXIT(1);
2803 }
2804
2805 outim = mri_new_vol_empty( ncol , nrow , 1 , MRI_double ) ;
2806 mri_fix_data_pointer( dtsar , outim ) ;
2807 mri_add_name( fname , outim ) ;
2808
2809 FRB(buf) ; RETURN(outim) ;
2810 }
2811
2812 /*---------------------------------------------------------------*/
2813
mri_read_complex_ascii(char * fname)2814 MRI_IMAGE * mri_read_complex_ascii( char * fname )
2815 {
2816 MRI_IMAGE * outim ;
2817 int ii,jj,val , used_tsar , alloc_tsar, ih ;
2818 float * tsar ;
2819 complex *ctsar;
2820 float temp ;
2821 FILE * fts ;
2822 char * ptr ;
2823 int ncol , bpos , blen , nrow ;
2824 static char *buf=NULL ; /* 20 Jun 2002: make a ptr */
2825
2826 floatvec *vec ; /* 20 Jul 2004 */
2827 int incts ;
2828
2829 ENTRY("mri_read_complex_ascii") ;
2830
2831 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
2832
2833 if( strncmp(fname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
2834 /*
2835 MRI_IMAGE *qim = mri_1D_complex_fromstring( fname+3 ) ;
2836 if( qim != NULL ){
2837 outim = mri_transpose(qim); mri_free(qim); RETURN(outim);
2838 }
2839 */
2840 fprintf(stderr,"Somebody was too lazy to allow this option here.\n"); RETURN(NULL);
2841 }
2842
2843 fts = fopen( fname , "r" );
2844 if( fts == NULL ){
2845 ERROR_message("mri_read_complex_ascii: couldn't open file %s" , fname ) ; RETURN(NULL) ;
2846 }
2847
2848 if( buf == NULL ) buf = AFMALL(char, LBUF) ; /* create buffer */
2849
2850 /** step 1: read in the first line and see how many numbers are in it
2851 (skipping lines that are comments or entirely blank) */
2852
2853 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2854 ptr = my_fgets( buf , LBUF , fts ) ;
2855 if( ptr==NULL || *ptr=='\0' ){ FRB(buf); fclose(fts); RETURN(NULL); } /* bad read? */
2856
2857 lbfill = 0.0f ; /* 10 Aug 2004 */
2858
2859 vec = decode_linebuf( buf ) ; /* 20 Jul 2004 */
2860 if( vec == NULL || vec->nar == 0 ){
2861 if (strlen(buf) > 64) {
2862 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2863 }
2864 if (linebufdied) {/* death? ZSS: Oct 19 2009*/
2865 static int nfail=0 ;
2866 if( ++nfail < 9 ) {
2867 fprintf(stderr,
2868 "\n** Error: Failed parsing data row 0 of 1D file '%.44s'\n"
2869 " Check for illegal non-numeric characters in:\n"
2870 " '%s'\n", fname,buf);
2871 if (nfail == 1)
2872 fprintf(stderr,
2873 " You can have text entries set to 0 with -ok_1D_text or by \n"
2874 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2875 }
2876 }
2877 if( vec != NULL ) KILL_floatvec(vec) ;
2878 FRB(buf); fclose(fts); RETURN(NULL);
2879 }
2880 ncol = vec->nar ; KILL_floatvec(vec) ;
2881 if (ncol % 2) {
2882 fprintf(stderr,"\n*** File does not have even number of columns."
2883 "\n That is a must for complex 1D files.\n");
2884 RETURN(NULL);
2885 }
2886 /* fprintf(stderr,"Have ncol = %d\n", ncol);*/
2887 /** At this point, ncol is the number of floats to be read from each line **/
2888
2889 rewind( fts ) ; /* will start over */
2890
2891 incts = MAX(INC_TSARSIZE,ncol) ;
2892 used_tsar = 0 ;
2893 alloc_tsar = incts ;
2894 tsar = (float *) malloc( sizeof(float) * alloc_tsar ) ;
2895 if( tsar == NULL ){
2896 fprintf(stderr,"\n*** malloc error in mri_read_complex_ascii ***\n"); EXIT(1);
2897 }
2898
2899 /** read lines, convert to floats, store **/
2900
2901 nrow = 0 ;
2902 while( 1 ){
2903 ptr = my_fgets( buf , LBUF , fts ) ; /* read */
2904 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
2905
2906 vec = decode_linebuf( buf ) ;
2907 if( vec == NULL ) break ;
2908 if( vec->nar == 0 ){ KILL_floatvec(vec); break; }
2909
2910 if( used_tsar + ncol >= alloc_tsar ){
2911 alloc_tsar += incts ;
2912 tsar = (float *)realloc( (void *)tsar,sizeof(float)*alloc_tsar );
2913 if( tsar == NULL ){
2914 fprintf(stderr,"\n*** realloc error in mri_read_complex_ascii ***\n"); EXIT(1);
2915 }
2916 }
2917 for( ii=0 ; ii < vec->nar && ii < ncol ; ii++ )
2918 tsar[used_tsar+ii] = vec->ar[ii] ;
2919 for( ; ii < ncol ; ii++ )
2920 tsar[used_tsar+ii] = 0.0 ;
2921 used_tsar += ncol ;
2922 KILL_floatvec(vec) ;
2923
2924 nrow++ ; /* got one more complete row! */
2925 }
2926 fclose( fts ) ; /* finished with this file! */
2927 (void) my_fgets( NULL , 0 , NULL ) ; /* reset [20 Jul 2004] */
2928
2929 if (linebufdied) {/* death? ZSS: Oct 19 2009 */
2930 static int nfail=0 ;
2931 if (strlen(buf) > 64) {
2932 buf[64]='\0'; buf[63]='.'; buf[62]='.'; buf[61]='.'; buf[60]=' ';
2933 }
2934 if( ++nfail < 9 ) {
2935 fprintf(stderr,
2936 "\n** Error: Failed parsing data row %d of 1D file '%.44s'\n"
2937 " Check for illegal non-numeric characters in:\n"
2938 " '%s'\n", nrow, fname,buf);
2939 if (nfail == 1)
2940 fprintf(stderr,
2941 " You can have text entries set to 0 with -ok_1D_text or by \n"
2942 " setting environment variable AFNI_1D_ZERO_TEXT to YES.\n");
2943 }
2944
2945 if (tsar) free(tsar); tsar = NULL;
2946 FRB(buf);
2947 RETURN(NULL);
2948 }
2949
2950 /* from <= 1 to < 1 (allow 1x1 image) 25 Jan 2006 [rickr] */
2951 if( used_tsar < 1 ){ FRB(buf); free(tsar); RETURN(NULL); }
2952
2953 tsar = (float *) realloc( tsar , sizeof(float) * used_tsar ) ;
2954 if( tsar == NULL ){
2955 fprintf(stderr,"\n*** final realloc error in mri_read_complex_ascii ***\n"); EXIT(1);
2956 }
2957
2958 /* now turn tsar into a complex vector */
2959 ctsar = (complex *) calloc(used_tsar, sizeof(complex));
2960 for( ii=0 ; ii < used_tsar; ii=ii+2) {
2961 /* fprintf(stderr,"tsar[%d]=%f\n", ii, tsar[ii]); */
2962 ih = ii/2;
2963 ctsar[ih].r = tsar[ii]; ctsar[ih].i = tsar[ii+1];
2964 }
2965
2966 outim = mri_new_vol_empty( ncol/2 , nrow , 1 , MRI_complex ) ;
2967 mri_fix_data_pointer( tsar , outim ) ;
2968 mri_add_name( fname , outim ) ;
2969
2970 FRB(buf) ; RETURN(outim) ;
2971 }
2972
2973 /*---------------------------------------------------------------------------*/
2974
2975 static char *dname=NULL ; static size_t ndname=0 ; /* 15 Nov 2007 */
2976
2977 #define DNAME_FIX(fn) \
2978 do{ size_t qq=strlen(fn)+7 ; \
2979 if( ndname < qq ){ dname=(char *)realloc((void *)dname,qq); ndname=qq; } \
2980 } while(0)
2981
2982 /*---------------------------------------------------------------------------*/
2983 /*! Read an ASCII file as columns, transpose to rows, allow column selectors.
2984
2985 \param fname = Input filename (max of 255 characters)
2986 \return Pointer to MRI_IMAGE if all went well; NULL if not.
2987 \date 16 Nov 1999
2988
2989 This function builds on mri_read_ascii() in two ways:
2990 - the input is transposed to rows (so that a 1x100 file becomes a 100x1 image)
2991 - column selectors [..] and row selectors {..} are allowed in fname
2992 - if fname ends in a ' character, file will be NOT be transposed
2993 */
2994
mri_read_1D(char * fname)2995 MRI_IMAGE * mri_read_1D( char *fname )
2996 {
2997 MRI_IMAGE *inim=NULL , *outim=NULL , *flim=NULL ;
2998 char *cpt , *dpt ;
2999 int ii,jj,nx,ny,nts , *ivlist , *ivl , *sslist ;
3000 float *far , *oar ;
3001 int flip ; /* 05 Sep 2006 */
3002
3003 ENTRY("mri_read_1D") ;
3004
3005 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
3006
3007 /*-- 14 Sep 2018: read a TSV file? --*/
3008
3009 cpt = strcasestr(fname,".tsv") ;
3010 if( cpt != NULL && ( cpt[4] == '\0' || cpt[4] == '[' ) ){
3011 NI_element *nel ;
3012 nel = THD_read_tsv(fname) ; /* cf. thd_table.c */
3013 if( nel != NULL ){
3014 outim = THD_niml_to_mri(nel) ; /* only numeric columns */
3015 NI_free_element(nel) ;
3016 if( outim != NULL ) RETURN(outim) ;
3017 } /* if it falls thru to here, read or conversion failed */
3018 }
3019
3020 /*-- 15 Apr 2019: read CSV? --*/
3021
3022 cpt = strcasestr(fname,".csv") ;
3023 if( cpt != NULL && ( cpt[4] == '\0' || cpt[4] == '[' ) ){
3024 NI_element *nel ;
3025 nel = THD_read_csv(fname) ; /* cf. thd_table.c */
3026 if( nel != NULL ){
3027 outim = THD_niml_to_mri(nel) ; /* only numeric columns */
3028 NI_free_element(nel) ;
3029 if( outim != NULL ) RETURN(outim) ;
3030 } /* if it falls thru to here, read or conversion failed */
3031 }
3032
3033 /*-- 25 Jan 2008: read from stdin? --*/
3034
3035 ii = strlen(fname) ;
3036 if( (ii <= 2 && fname[0] == '-') ||
3037 (ii <= 6 && strncmp(fname,"stdin" ,5) == 0) ||
3038 (ii <= 9 && strncmp(fname,"1D:stdin",8) == 0) || /* Isaac S.'s 02/06/13 */
3039 (ii <= 9 && strncmp(fname,"/dev/fd0",8) == 0) ){
3040 inim = mri_read_1D_stdin() ;
3041 if( inim != NULL && fname[ii-1] == '\'' ){
3042 flim = mri_transpose(inim); mri_free(inim); inim = flim;
3043 } else if( inim == NULL ){
3044 ERROR_message("mri_read_1D: can't read 1D data from stdin") ;
3045 }
3046 RETURN(inim) ;
3047 }
3048
3049 if( strncasecmp(fname,"jRandom1D:",10) == 0 && isdigit(fname[10]) ){ /* 17 Mar 2016 */
3050 int nx=0 , ny=0 ;
3051 if( strchr(fname+10,':') != NULL )
3052 sscanf( fname+10 , "%d:%d" , &nx,&ny ) ;
3053 else
3054 sscanf( fname+10 , "%d,%d" , &nx,&ny ) ;
3055 inim = jRandom1D(nx,ny) ;
3056 if( inim == NULL )
3057 WARNING_message("Can't decode %s",fname) ;
3058 RETURN(inim) ;
3059 }
3060
3061 /* read from a 3D: file? [31 Aug 2018] */
3062
3063 if( strlen(fname) > 9 &&
3064 fname[0] == '3' &&
3065 fname[1] == 'D' &&
3066 (fname[2] == ':' || fname[3] == ':') ){
3067
3068 MRI_IMARR *imar = mri_read_3D(fname) ;
3069 if( imar == NULL ){
3070 ERROR_message("mri_read_1D: can't read valid data from %s",fname) ;
3071 RETURN(NULL);
3072 }
3073 if( IMARR_COUNT(imar) == 0 ){ DESTROY_IMARR(imar); RETURN(NULL); }
3074 outim = mri_to_float( IMARR_SUBIM(imar,0) ) ;
3075 DESTROY_IMARR(imar) ; RETURN(outim) ;
3076 }
3077
3078 /*-- back to reading from an actual file --*/
3079
3080 DNAME_FIX(fname) ;
3081 strcpy(dname,fname); ii = strlen(dname); /* 05 Sep 2006 */
3082 flip = (dname[ii-1] == '\''); if( flip ) dname[ii-1] = '\0';
3083
3084 if( strncmp(dname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
3085 outim = mri_1D_fromstring( dname+3 ) ;
3086 if( outim == NULL ){
3087 ERROR_message("mri_read_1D: can't read from string '%s'",dname) ;
3088 } else {
3089 if( flip ){ inim=mri_transpose(outim); mri_free(outim); outim=inim; }
3090 mri_add_name("1D:...",outim) ;
3091 }
3092 RETURN(outim) ;
3093 }
3094
3095 /*-- split filename and subvector list --*/
3096
3097 cpt = strstr(fname,"[") ;
3098 dpt = strstr(fname,"{") ; /* 30 Apr 2003: subsampling list */
3099
3100 if( cpt == fname || dpt == fname ){ /* can't be at start of filename! */
3101 ERROR_message("Illegal filename in mri_read_1D('%s')\n",fname) ;
3102 RETURN(NULL) ;
3103 } else { /* got a subvector list */
3104 if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; }
3105 if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; }
3106 ii = strlen(dname) ;
3107 if( dname[ii-1] == '\'' ){ flip = 0 ; dname[ii-1] = '\0' ; } /* 16 Jun 2010 */
3108 }
3109
3110 /*-- read file in, flip it sideways --*/
3111
3112 if( strcmp(dname,"stdin") != 0 && strncmp(dname,"/dev/fd0",8) != 0 ){
3113
3114 { FILE *qp = fopen(dname,"r") ; /* from a named pipe/fifo [26 Aug 2018] */
3115 if( qp != NULL ){
3116 struct stat qst ; int qq ;
3117 qq = fstat( fileno(qp) , &qst ) ;
3118 if( qq == 0 && S_ISFIFO(qst.st_mode) ){ /* it IS a fifo! */
3119 flim = mri_read_1D_pipe(qp) ; /* don't flip it */
3120 }
3121 fclose(qp) ; /* done with this, no matter what happened */
3122 }
3123 }
3124 if( flim == NULL ){ /* try reading as normal file */
3125 inim = mri_read_ascii(dname) ;
3126 if( inim == NULL ) RETURN(NULL) ;
3127 flim = mri_transpose(inim) ; mri_free(inim) ;
3128 }
3129 } else {
3130 flim = mri_read_1D_stdin() ; /* stdin with subvector selection */
3131 if( flim == NULL ) RETURN(NULL);
3132 }
3133
3134 /*-- get the subvector and subsampling lists, if any --*/
3135
3136 nx = flim->nx ; ny = flim->ny ;
3137
3138 ivlist = MCW_get_intlist( ny , cpt ) ; /* subvector list */
3139 sslist = MCW_get_intlist( nx , dpt ) ; /* subsampling list */
3140
3141 /* if either columns(vectors) or rows (subsamples) are specified */
3142 /* make sure the lists are valid. Return null if not */
3143 if( ((cpt!= NULL) && (ivlist==NULL)) ||
3144 ((dpt!= NULL) && (sslist==NULL)) )
3145 RETURN(NULL);
3146
3147 /* if have subvector list, extract those rows into a new image */
3148
3149 if( ivlist != NULL && ivlist[0] > 0 ){
3150 nts = ivlist[0] ; /* number of subvectors */
3151 ivl = ivlist + 1 ; /* start of array of subvectors */
3152
3153 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3154 if( ivl[ii] < 0 || ivl[ii] >= ny ){
3155 fprintf(stderr,"*** Out-of-range subvector [list] in mri_read_1D: %s\n",fname) ;
3156 mri_free(flim) ; free(ivlist) ; RETURN(NULL) ;
3157 }
3158 }
3159
3160 outim = mri_new( nx , nts , MRI_float ) ; /* make output image */
3161 far = MRI_FLOAT_PTR( flim ) ;
3162 oar = MRI_FLOAT_PTR( outim ) ;
3163
3164 for( ii=0 ; ii < nts ; ii++ ) /* copy desired rows */
3165 memcpy( oar + ii*nx , far + ivl[ii]*nx , sizeof(float)*nx ) ;
3166
3167 mri_free(flim); free(ivlist); flim = outim; ny = nts;
3168 }
3169
3170 /* if have subsampling list, extract those columns into a new image */
3171
3172 if( sslist != NULL && sslist[0] > 0 ){
3173 nts = sslist[0] ; /* number of columns to get */
3174 ivl = sslist + 1 ; /* start of array of column indexes */
3175
3176 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3177 if( ivl[ii] < 0 || ivl[ii] >= nx ){
3178 fprintf(stderr,"*** Out-of-range subsampling {list} in mri_read_1D: %s\n",fname) ;
3179 mri_free(flim) ; free(sslist) ; RETURN(NULL) ;
3180 }
3181 }
3182
3183 if( AFNI_yesenv("AFNI_TEST_SUBSETTER") ){
3184 INFO_message("Calling mri_subset_x2D") ;
3185 outim = mri_subset_x2D( nts , ivl , flim ) ;
3186 } else {
3187 outim = mri_new( nts , ny , MRI_float ) ; /* make output image */
3188 far = MRI_FLOAT_PTR( flim ) ;
3189 oar = MRI_FLOAT_PTR( outim ) ;
3190
3191 for( ii=0 ; ii < nts ; ii++ ) /* copy desired columns */
3192 for( jj=0 ; jj < ny ; jj++ )
3193 oar[ii+jj*nts] = far[ivl[ii]+jj*nx] ;
3194 }
3195
3196 mri_free(flim); free(sslist); flim = outim;
3197 }
3198
3199 if( flip ){ inim=mri_transpose(flim); mri_free(flim); flim=inim; }
3200
3201 mri_add_name(fname,flim) ; RETURN(flim) ;
3202 }
3203
3204 /*---------------------------------------------------------------------------*/
3205
mri_rowmajorize_1D(MRI_IMAGE * im)3206 MRI_IMAGE * mri_rowmajorize_1D( MRI_IMAGE *im )
3207 {
3208 MRI_IMAGE *qim = mri_transpose(im) ;
3209 qim->nx *= qim->ny ;
3210 qim->ny = 1 ;
3211 return qim ;
3212 }
3213
3214 /*---------------------------------------------------------------------------*/
3215 /*! Read header lines from a 1D file into a newly-malloc()-ed string. */
3216
mri_read_1D_headerlines(char * fname)3217 char * mri_read_1D_headerlines( char *fname )
3218 {
3219 int ii , nout=0 ;
3220 FILE *fp ;
3221 char lbuf[NLL] , *cout=NULL , *dpt ;
3222
3223 ENTRY("mri_read_1D_headerlines") ;
3224
3225 if( fname == NULL || *fname == '\0' ) RETURN(NULL) ;
3226 if( strncmp(fname,"1D:",3) == 0 ) RETURN(NULL) ;
3227
3228 ii = strlen(fname) ;
3229 if( (ii <= 2 && fname[0] == '-') ||
3230 (ii <= 6 && strncmp(fname,"stdin" ,5) == 0) ||
3231 (ii <= 9 && strncmp(fname,"/dev/fd0",8) == 0) ){
3232
3233 fp = stdin ;
3234 } else {
3235 fp = fopen( fname , "r" ) ;
3236 if( fp == NULL ){
3237 ERROR_message("mri_read_1D_headerlines: couldn't open file %s" , fname ) ; RETURN(NULL) ;
3238 }
3239 }
3240
3241 /* read # lines, catenate them */
3242
3243 while(1){
3244 lbuf[0] = '\0' ;
3245 dpt = afni_fgets( lbuf , NLL , fp ) ; /* read a line of data */
3246 if( dpt == NULL ) break ; /* nothing more to read */
3247 ii = strlen(lbuf) ; if( ii == 0 ) break ; /* nada => finished */
3248 if( lbuf[0] != '#' ) break ; /* not '#' => finished */
3249 cout = (char *)realloc( cout , sizeof(char)*(nout+ii+2) ) ; /* space */
3250 strcpy( cout+nout , lbuf ) ; /* copy new stuff into space */
3251 nout = strlen(cout) ; /* length of output now */
3252 }
3253 if( fp != stdin ) fclose(fp) ;
3254
3255 RETURN(cout) ;
3256 }
3257
3258 /*---------------------------------------------------------------------------*/
3259 /*! Read a 1D with a total of 12 or 16 numbers and return it as a 4x4 matrix
3260 for spatial affine transformation
3261
3262 \param fname = Input filename (max of 255 characters)
3263 \return Pointer to MRI_IMAGE if all went well; NULL if not.
3264 \date Nov 24 2009
3265
3266 See mri_read_1D for special name modifiers.
3267
3268 The function accepts:
3269 one row, or one column vector of 12 or 16 elements
3270 one 3x4, or one 4x4 matrix
3271
3272 When a total of 16 elements are read in, the last row
3273 must be 0 0 0 1.
3274
3275 Transposing is allowed but I don't know why you'd want it.
3276
3277 stdin input is allowed
3278 */
3279
mri_read_4x4AffXfrm_1D(char * fname)3280 MRI_IMAGE * mri_read_4x4AffXfrm_1D( char *fname )
3281 {
3282 MRI_IMAGE *inim =NULL , *outim =NULL;
3283 char *cpt , *dpt ;
3284 int ii=0, c=0, r=0 ;
3285 float *far , *oar ;
3286 int flip=0; /* 05 Sep 2006 */
3287
3288 ENTRY("mri_read_4x4AffXfrm_1D") ;
3289
3290 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
3291
3292 /*-- 25 Jan 2008: read from stdin? --*/
3293
3294 ii = strlen(fname) ;
3295 if( (ii <= 2 && fname[0] == '-') ||
3296 (ii <= 6 && strncmp(fname,"stdin" ,5) == 0) ||
3297 (ii <= 9 && strncmp(fname,"/dev/fd0",8) == 0) ){
3298 inim = mri_read_1D_stdin() ;
3299 if (!inim) RETURN(inim);
3300 if( inim != NULL && fname[ii-1] == '\'' ){
3301 flip = 1; /*transpose later */
3302 }
3303 } else {
3304 /*-- Read from file, transposing is handled later --*/
3305 DNAME_FIX(fname) ;
3306 strcpy(dname,fname); ii = strlen(dname); /* 05 Sep 2006 */
3307 flip = (dname[ii-1] == '\''); if( flip ) dname[ii-1] = '\0';
3308
3309 /* read with regular function, no transposing yet */
3310 inim = mri_read_1D(dname);
3311 if (!inim) RETURN(inim);
3312 }
3313
3314 /* check on dimensions */
3315 if ( !(inim->nx == 1 && inim->ny == 12) &&
3316 !(inim->nx == 12 && inim->ny == 1 ) &&
3317 !(inim->nx == 1 && inim->ny == 16) &&
3318 !(inim->nx == 16 && inim->ny == 1 ) &&
3319 !(inim->nx == 3 && inim->ny == 4 ) &&
3320 !(inim->nx == 4 && inim->ny == 4 ) ) {
3321 fprintf( stderr,
3322 "*** Bad dimensions of %dx%d in mri_read_4x4AffXfmr_1D: %s\n"
3323 " Allowed dimensions are 12x1, 1x12, 16x1, 1x16, 3x4, and 4x4\n",
3324 inim->nx, inim->ny, dname) ;
3325 mri_free(inim); inim = NULL;
3326 RETURN(inim);
3327 }
3328
3329 /* fprintf(stderr,"%d x %d\n", inim->nx, inim->ny); */
3330
3331 /* prepare output */
3332 outim = mri_new( 4 , 4 , MRI_float ) ; /* make output image */
3333 far = MRI_FLOAT_PTR( inim ) ;
3334 oar = MRI_FLOAT_PTR( outim ) ;
3335
3336 /* fillup oar, column by column */
3337 if ((inim->nx == 1 && inim->ny == 12) ||
3338 (inim->nx == 12 && inim->ny == 1 ) ) {
3339 ii = 0;
3340 for (r=0; r<4; ++r) {
3341 for (c=0; c<4; ++c) {
3342 if (r < 3) {
3343 oar[r+c*4] = far[ii]; ++ii;
3344 } else {
3345 oar[r+c*4] = 0.0;
3346 }
3347 }
3348 }
3349 oar[15] = 1.0;
3350 } else if ((inim->nx == 3 && inim->ny == 4 ) ) {
3351 for (r=0; r<4; ++r) {
3352 for (c=0; c<4; ++c) {
3353 if (r < 3) {
3354 oar[r+c*4] = far[r+c*3];
3355 } else {
3356 oar[r+c*4] = 0.0;
3357 }
3358 }
3359 }
3360 oar[15] = 1.0;
3361 } else if ( (inim->nx == 1 && inim->ny == 16) ||
3362 (inim->nx == 16 && inim->ny == 1 ) ) {
3363 ii = 0;
3364 for (r=0; r<4; ++r) {
3365 for (c=0; c<4; ++c) {
3366 oar[r+c*4] = far[ii]; ++ii;
3367 }
3368 }
3369 } else if ( (inim->nx == 4 && inim->ny == 4 ) ) {
3370 for (r=0; r<4; ++r) {
3371 for (c=0; c<4; ++c) {
3372 oar[r+c*4] = far[r+c*4];
3373 }
3374 }
3375 }
3376
3377 /* all done with input */
3378 mri_free(inim); inim = NULL;
3379
3380 /* final check */
3381 if ( oar[15] != 1.0f ||
3382 oar[11] != 0.0f ||
3383 oar[ 7] != 0.0f ||
3384 oar[ 3] != 0.0f ) {
3385 fprintf( stderr,
3386 "*** Bad 4th row values in %s.\n"
3387 "Expecting 0.0 0.0 0.0 1.0, \n"
3388 "got %f %f %f %f\n",
3389 dname, oar[3], oar[7], oar[11], oar[15]);
3390 mri_free(outim); outim = NULL;
3391 RETURN(outim);
3392 }
3393
3394 if( flip ){
3395 fprintf(stderr,"Transposing xform!\n");
3396 inim=mri_transpose(outim);
3397 mri_free(outim); outim=inim; inim=NULL;
3398 }
3399
3400 mri_add_name(fname,outim) ; RETURN(outim) ;
3401 }
3402
mri_read_double_1D(char * fname)3403 MRI_IMAGE * mri_read_double_1D( char *fname )
3404 {
3405 MRI_IMAGE *inim , *outim , *flim ;
3406 char *cpt , *dpt ;
3407 int ii,jj,nx,ny,nts , *ivlist , *ivl , *sslist ;
3408 double *dar , *oar ;
3409 int flip ; /* 05 Sep 2006 */
3410
3411 ENTRY("mri_read_double_1D") ;
3412
3413 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
3414 DNAME_FIX(fname) ;
3415 strcpy(dname,fname); ii = strlen(dname); /* 05 Sep 2006 */
3416 flip = (dname[ii-1] == '\''); if( flip ) dname[ii-1] = '\0';
3417
3418 if( strncmp(dname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
3419 /*
3420 outim = mri_1D_double_fromstring( dname+3 ) ;
3421 if( flip ){ inim=mri_transpose(outim); mri_free(outim); outim=inim; }
3422 RETURN(outim) ;
3423 */
3424 ERROR_message("Somebody nameless was too lazy to allow this option here."); RETURN(NULL);
3425 }
3426
3427 /*-- split filename and subvector list --*/
3428
3429 cpt = strstr(fname,"[") ;
3430 dpt = strstr(fname,"{") ; /* 30 Apr 2003: subsampling list */
3431
3432 if( cpt == fname || dpt == fname ){ /* can't be at start of filename! */
3433 ERROR_message("Illegal filename in mri_read_double_1D('%s')\n",fname) ;
3434 RETURN(NULL) ;
3435 } else { /* got a subvector list */
3436 if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; }
3437 if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; }
3438 }
3439
3440 /*-- read file in, flip it sideways --*/
3441
3442 inim = mri_read_double_ascii(dname) ;
3443 if( inim == NULL ) RETURN(NULL) ;
3444 flim = mri_transpose(inim) ; mri_free(inim) ;
3445 if( flim == NULL ) {
3446 fprintf(stderr, "Failed to transpose image\n");
3447 RETURN(NULL) ;
3448 }
3449 /*-- get the subvector and subsampling lists, if any --*/
3450
3451 nx = flim->nx ; ny = flim->ny ;
3452
3453 ivlist = MCW_get_intlist( ny , cpt ) ; /* subvector list */
3454 sslist = MCW_get_intlist( nx , dpt ) ; /* subsampling list */
3455
3456 /* if have subvector list, extract those rows into a new image */
3457
3458 if( ivlist != NULL && ivlist[0] > 0 ){
3459 nts = ivlist[0] ; /* number of subvectors */
3460 ivl = ivlist + 1 ; /* start of array of subvectors */
3461
3462 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3463 if( ivl[ii] < 0 || ivl[ii] >= ny ){
3464 fprintf(stderr,"*** Out-of-range subvector [list] in mri_read_double_1D: %s\n",fname) ;
3465 mri_free(flim) ; free(ivlist) ; RETURN(NULL) ;
3466 }
3467 }
3468
3469 outim = mri_new( nx , nts , MRI_double ) ; /* make output image */
3470 dar = MRI_DOUBLE_PTR( flim ) ;
3471 oar = MRI_DOUBLE_PTR( outim ) ;
3472
3473 for( ii=0 ; ii < nts ; ii++ ) /* copy desired rows */
3474 memcpy( oar + ii*nx , dar + ivl[ii]*nx , sizeof(double)*nx ) ;
3475
3476 mri_free(flim); free(ivlist); flim = outim; ny = nts;
3477 }
3478
3479 /* if have subsampling list, extract those columns into a new image */
3480
3481 if( sslist != NULL && sslist[0] > 0 ){
3482 nts = sslist[0] ; /* number of columns to get */
3483 ivl = sslist + 1 ; /* start of array of column indexes */
3484
3485 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3486 if( ivl[ii] < 0 || ivl[ii] >= nx ){
3487 fprintf(stderr,"*** Out-of-range subsampling {list} in mri_read_double_1D: %s\n",fname) ;
3488 mri_free(flim) ; free(sslist) ; RETURN(NULL) ;
3489 }
3490 }
3491
3492 outim = mri_new( nts , ny , MRI_double ) ; /* make output image */
3493 dar = MRI_DOUBLE_PTR( flim ) ;
3494 oar = MRI_DOUBLE_PTR( outim ) ;
3495
3496 for( ii=0 ; ii < nts ; ii++ ) /* copy desired columns */
3497 for( jj=0 ; jj < ny ; jj++ )
3498 oar[ii+jj*nts] = dar[ivl[ii]+jj*nx] ;
3499
3500 mri_free(flim); free(sslist); flim = outim;
3501 }
3502
3503 if( flip ){ inim=mri_transpose(flim); mri_free(flim); flim=inim; }
3504
3505 mri_add_name(fname,flim) ; RETURN(flim) ;
3506 }
3507
mri_read_complex_1D(char * fname)3508 MRI_IMAGE * mri_read_complex_1D( char *fname )
3509 {
3510 MRI_IMAGE *inim , *outim , *flim ;
3511 char *cpt , *dpt ;
3512 int ii,jj,nx,ny,nts , *ivlist , *ivl , *sslist ;
3513 complex *far , *oar ;
3514 int flip ; /* 05 Sep 2006 */
3515
3516 ENTRY("mri_read_complex_1D") ;
3517
3518 if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
3519 DNAME_FIX(fname) ;
3520 strcpy(dname,fname); ii = strlen(dname); /* 05 Sep 2006 */
3521 flip = (dname[ii-1] == '\''); if( flip ) dname[ii-1] = '\0';
3522
3523 if( strncmp(dname,"1D:",3) == 0 ){ /* 28 Apr 2003 */
3524 /*
3525 outim = mri_1D_complex_fromstring( dname+3 ) ;
3526 if( flip ){ inim=mri_transpose(outim); mri_free(outim); outim=inim; }
3527 RETURN(outim) ;
3528 */
3529 ERROR_message("Somebody was too lazy to allow this option here."); RETURN(NULL);
3530 }
3531
3532 /*-- split filename and subvector list --*/
3533
3534 cpt = strstr(fname,"[") ;
3535 dpt = strstr(fname,"{") ; /* 30 Apr 2003: subsampling list */
3536
3537 if( cpt == fname || dpt == fname ){ /* can't be at start of filename! */
3538 ERROR_message("Illegal filename in mri_read_complex_1D('%s')\n",fname) ;
3539 RETURN(NULL) ;
3540 } else { /* got a subvector list */
3541 if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; }
3542 if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; }
3543 }
3544
3545 /*-- read file in, flip it sideways --*/
3546
3547 inim = mri_read_complex_ascii(dname) ;
3548 if( inim == NULL ) RETURN(NULL) ;
3549 flim = mri_transpose(inim) ; mri_free(inim) ;
3550
3551 /*-- get the subvector and subsampling lists, if any --*/
3552
3553 nx = flim->nx ; ny = flim->ny ;
3554
3555 ivlist = MCW_get_intlist( ny , cpt ) ; /* subvector list */
3556 sslist = MCW_get_intlist( nx , dpt ) ; /* subsampling list */
3557
3558 /* if have subvector list, extract those rows into a new image */
3559
3560 if( ivlist != NULL && ivlist[0] > 0 ){
3561 nts = ivlist[0] ; /* number of subvectors */
3562 ivl = ivlist + 1 ; /* start of array of subvectors */
3563
3564 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3565 if( ivl[ii] < 0 || ivl[ii] >= ny ){
3566 fprintf(stderr,"*** Out-of-range subvector [list] in mri_read_complex_1D: %s\n",fname) ;
3567 mri_free(flim) ; free(ivlist) ; RETURN(NULL) ;
3568 }
3569 }
3570
3571 outim = mri_new( nx , nts , MRI_complex ) ; /* make output image */
3572 far = MRI_COMPLEX_PTR( flim ) ;
3573 oar = MRI_COMPLEX_PTR( outim ) ;
3574
3575 for( ii=0 ; ii < nts ; ii++ ) /* copy desired rows */
3576 memcpy( oar + ii*nx , far + ivl[ii]*nx , sizeof(complex)*nx ) ;
3577
3578 mri_free(flim); free(ivlist); flim = outim; ny = nts;
3579 }
3580
3581 /* if have subsampling list, extract those columns into a new image */
3582
3583 if( sslist != NULL && sslist[0] > 0 ){
3584 nts = sslist[0] ; /* number of columns to get */
3585 ivl = sslist + 1 ; /* start of array of column indexes */
3586
3587 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */
3588 if( ivl[ii] < 0 || ivl[ii] >= nx ){
3589 fprintf(stderr,"*** Out-of-range subsampling {list} in mri_read_complex_1D: %s\n",fname) ;
3590 mri_free(flim) ; free(sslist) ; RETURN(NULL) ;
3591 }
3592 }
3593
3594 outim = mri_new( nts , ny , MRI_complex ) ; /* make output image */
3595 far = MRI_COMPLEX_PTR( flim ) ;
3596 oar = MRI_COMPLEX_PTR( outim ) ;
3597
3598 for( ii=0 ; ii < nts ; ii++ ) /* copy desired columns */
3599 for( jj=0 ; jj < ny ; jj++ ) {
3600 oar[ii+jj*nts].r = far[ivl[ii]+jj*nx].r ;
3601 oar[ii+jj*nts].i = far[ivl[ii]+jj*nx].i ;
3602 }
3603
3604 mri_free(flim); free(sslist); flim = outim;
3605 }
3606
3607 if( flip ){ inim=mri_transpose(flim); mri_free(flim); flim=inim; }
3608
3609 mri_add_name(fname,flim) ; RETURN(flim) ;
3610 }
3611
3612 /*-----------------------------------------------------------------------------------*/
3613
3614 static MRI_IMAGE *im_stdin = NULL ;
3615
mri_copy_1D_stdin(void)3616 MRI_IMAGE * mri_copy_1D_stdin(void) /* 05 Mar 2010 */
3617 {
3618 if( im_stdin == NULL ) im_stdin = mri_read_1D_stdin() ;
3619 return mri_copy(im_stdin) ;
3620 }
3621
mri_clear_1D_stdin(void)3622 void mri_clear_1D_stdin(void)
3623 {
3624 if( im_stdin != NULL ){ mri_free(im_stdin); im_stdin = NULL; }
3625 }
3626
3627 /*-----------------------------------------------------------------------------------*/
3628 /* Read a 1D file from stdin; adapted from 1dplot.c */
3629
mri_read_1D_stdin(void)3630 MRI_IMAGE * mri_read_1D_stdin(void)
3631 {
3632 #if 1
3633 return ( mri_read_1D_pipe(stdin) ) ; /* the new way [26 Aug 2019] */
3634
3635 #else /*---------- the olden way (stdin specific) ----------*/
3636 #define SIN_NLBUF 131072
3637 #define SIN_NVMAX 10000
3638 char *lbuf , *cpt , *dpt ;
3639 int nval , ii,nx,ny ;
3640 float *val , fff , *far ;
3641 MRI_IMAGE *flim , *inim ;
3642
3643 ENTRY("mri_read_1D_stdin") ;
3644
3645 if( im_stdin != NULL ){
3646 ININFO_message("copying im_stdin") ;
3647 inim = mri_copy(im_stdin); RETURN(inim);
3648 }
3649
3650 /*** INFO_message("reading 1D_stdin") ; ***/
3651 lbuf = (char * )malloc(sizeof(char )*SIN_NLBUF) ;
3652 val = (float *)malloc(sizeof(float)*SIN_NVMAX) ;
3653
3654 do{ /* read lines until 1st char is non-blank and non-# */
3655 cpt = afni_fgets(lbuf,SIN_NLBUF,stdin) ;
3656 if( cpt==NULL ){ free(val);free(lbuf); RETURN(NULL); }
3657 for( ii=0 ; cpt[ii] != '\0' && isspace(cpt[ii]) ; ii++ ) ; /* nada */
3658 } while( cpt[ii] == '\0' || cpt[ii] == '#' ) ;
3659
3660 nval = 0 ; cpt = lbuf ; /* read numbers from lbuf into val */
3661 while(1){
3662 fff = strtod(cpt,&dpt) ; if( dpt == cpt ) break ;
3663 val[nval++] = fff ; if( nval == SIN_NVMAX ) break ;
3664 cpt = dpt; if( *cpt == ',' ) cpt++; if( *cpt == '\0' ) break;
3665 }
3666 if( nval < 1 ){ free(val);free(lbuf); RETURN(NULL); }
3667
3668 nx = nval ; ny = 1 ;
3669 far = (float *) malloc(sizeof(float)*nx) ;
3670 memcpy(far,val,sizeof(float)*nx) ;
3671
3672 while(1){ /* read from stdin */
3673 cpt = afni_fgets(lbuf,SIN_NLBUF,stdin) ;
3674 if( cpt == NULL ) break ; /* done */
3675 for( ii=0 ; cpt[ii] != '\0' && isspace(cpt[ii]) ; ii++ ) ; /* nada */
3676 if( cpt[ii] == '\0' || cpt[ii] == '#' ) continue ; /* skip */
3677
3678 memset(val,0,sizeof(float)*nx) ; /* set input buffer to zero */
3679 nval = 0 ; cpt = lbuf ; /* read numbers from lbuf into val */
3680 while(1){
3681 fff = strtod(cpt,&dpt) ; if( dpt == cpt ) break ;
3682 val[nval++] = fff ; if( nval == nx ) break ;
3683 cpt = dpt; if( *cpt == ',' ) cpt++; if( *cpt == '\0' ) break;
3684 }
3685 far = (float *) realloc( far , sizeof(float)*(ny+1)*nx ) ;
3686 memcpy(far+ny*nx,val,sizeof(float)*nx) ; ny++ ;
3687 }
3688
3689 flim = mri_new_vol_empty( nx,ny,1 , MRI_float ) ;
3690 mri_fix_data_pointer( far , flim ) ;
3691 if( ny > 1 ){ /* more than one row ==> transpose (the usual case) */
3692 inim = mri_transpose(flim) ; mri_free(flim) ;
3693 } else { /* only 1 row ==> am OK this way */
3694 inim = flim ;
3695 }
3696 free((void *)val); free((void *)lbuf);
3697 mri_add_name("stdin",inim); im_stdin = mri_copy(inim); RETURN(inim);
3698 #endif /*---------- end of the olden way ----------*/
3699
3700 }
3701
3702 /*-----------------------------------------------------------------------------------*/
3703 /* Read a 1D file from a stream with no seek/rewind (e.g., a pipe) */
3704
mri_read_1D_pipe(FILE * fp)3705 MRI_IMAGE * mri_read_1D_pipe( FILE *fp )
3706 {
3707 #define SIN_NLBUF 131072
3708 #define SIN_NVMAX 10000
3709 char *lbuf , *cpt , *dpt ;
3710 int nval , ii,nx,ny ;
3711 float *val , fff , *far ;
3712 MRI_IMAGE *flim , *inim ;
3713
3714 ENTRY("mri_read_1D_pipe") ;
3715
3716 if( fp == NULL ) RETURN(NULL) ;
3717
3718 if( fp == stdin && im_stdin != NULL ){
3719 ININFO_message("copying im_stdin") ;
3720 inim = mri_copy(im_stdin); RETURN(inim);
3721 }
3722
3723 lbuf = (char * )malloc(sizeof(char )*SIN_NLBUF) ;
3724 val = (float *)malloc(sizeof(float)*SIN_NVMAX) ;
3725
3726 do{ /* read lines until 1st char is non-blank and non-# */
3727 cpt = afni_fgets(lbuf,SIN_NLBUF,fp) ;
3728 if( cpt==NULL ){ free(val);free(lbuf); RETURN(NULL); }
3729 for( ii=0 ; cpt[ii] != '\0' && isspace(cpt[ii]) ; ii++ ) ; /* nada */
3730 } while( cpt[ii] == '\0' || cpt[ii] == '#' ) ;
3731
3732 nval = 0 ; cpt = lbuf ; /* read numbers from lbuf into val */
3733 while(1){
3734 fff = strtod(cpt,&dpt) ; if( dpt == cpt ) break ;
3735 val[nval++] = fff ; if( nval == SIN_NVMAX ) break ;
3736 cpt = dpt; if( *cpt == ',' ) cpt++; if( *cpt == '\0' ) break;
3737 }
3738 if( nval < 1 ){ free(val);free(lbuf); RETURN(NULL); }
3739
3740 nx = nval ; ny = 1 ;
3741 far = (float *) malloc(sizeof(float)*nx) ;
3742 memcpy(far,val,sizeof(float)*nx) ;
3743
3744 while(1){ /* read from fp */
3745 cpt = afni_fgets(lbuf,SIN_NLBUF,fp) ;
3746 if( cpt == NULL ) break ; /* done */
3747 for( ii=0 ; cpt[ii] != '\0' && isspace(cpt[ii]) ; ii++ ) ; /* nada */
3748 if( cpt[ii] == '\0' || cpt[ii] == '#' ) continue ; /* skip */
3749
3750 memset(val,0,sizeof(float)*nx) ; /* set input buffer to zero */
3751 nval = 0 ; cpt = lbuf ; /* read numbers from lbuf into val */
3752 while(1){
3753 fff = strtod(cpt,&dpt) ; if( dpt == cpt ) break ;
3754 val[nval++] = fff ; if( nval == nx ) break ;
3755 cpt = dpt; if( *cpt == ',' ) cpt++; if( *cpt == '\0' ) break;
3756 }
3757 far = (float *) realloc( far , sizeof(float)*(ny+1)*nx ) ;
3758 memcpy(far+ny*nx,val,sizeof(float)*nx) ; ny++ ;
3759 }
3760
3761 flim = mri_new_vol_empty( nx,ny,1 , MRI_float ) ;
3762 mri_fix_data_pointer( far , flim ) ;
3763 if( ny > 1 ){ /* more than one row ==> transpose (the usual case) */
3764 inim = mri_transpose(flim) ; mri_free(flim) ;
3765 } else { /* only 1 row ==> am OK this way */
3766 inim = flim ;
3767 }
3768 free((void *)val); free((void *)lbuf);
3769 if( fp == stdin ){
3770 if( im_stdin == NULL ) mri_free(im_stdin) ;
3771 mri_add_name("stdin",inim); im_stdin = mri_copy(inim);
3772 }
3773 RETURN(inim);
3774 }
3775
3776 /*-----------------------------------------------------------------------------------*/
3777 /* Read ragged rows, with '*' being set to the filler value [28 Jul 2004] */
3778
mri_read_ascii_ragged(char * fname,float filler)3779 MRI_IMAGE * mri_read_ascii_ragged( char *fname , float filler )
3780 {
3781 MRI_IMAGE *outim ;
3782 int ii,jj , ncol,nrow ;
3783 float *tsar ;
3784 FILE *fts ;
3785 char *ptr ;
3786 static char *buf=NULL ;
3787 floatvec *fvec ;
3788
3789 ENTRY("mri_read_ascii_ragged") ;
3790
3791 if( fname == NULL || *fname == '\0' ){ FRB(buf); RETURN(NULL); }
3792
3793 if( strncmp(fname,"1D:",3) == 0 ){ /* 05 Jan 2007 */
3794 outim = mri_read_ragged_fromstring( fname+3 , filler ) ;
3795 FRB(buf); RETURN(outim) ;
3796 }
3797
3798 fts = fopen( fname , "r" );
3799 if( fts == NULL ){
3800 FRB(buf) ;
3801 ERROR_message("mri_read_ascii_ragged: couldn't open file %s" , fname ) ; RETURN(NULL) ;
3802 }
3803
3804 if( buf == NULL ) buf = AFMALL(char, LBUF) ;
3805
3806 /** step 1: read in ALL lines, see how many numbers are in each,
3807 in order to get the maximum row length and # of rows **/
3808
3809 lbfill = filler ; /* 10 Aug 2004 */
3810
3811 (void) my_fgets( NULL , 0 , NULL ) ; /* reset */
3812 ncol = nrow = 0 ;
3813 while(1){
3814 ptr = my_fgets( buf , LBUF , fts ) ;
3815 if( ptr==NULL || *ptr=='\0' ) break ;
3816 fvec = decode_linebuf( buf ) ;
3817 if( fvec != NULL && fvec->nar > 0 ){ nrow++; ncol = MAX(ncol,fvec->nar); }
3818 if( fvec != NULL ) KILL_floatvec(fvec) ; else break ;
3819 }
3820 if( nrow == 0 || ncol == 0 ){ fclose(fts); FRB(buf); lbfill=0.0f; RETURN(NULL); }
3821
3822 /** At this point, ncol is the number of floats to be read from each line **/
3823
3824 rewind( fts ) ; /* will start over */
3825
3826 outim = mri_new( ncol , nrow , MRI_float ) ;
3827 tsar = MRI_FLOAT_PTR(outim) ;
3828
3829 /** read lines, convert to floats, store **/
3830
3831 nrow = 0 ;
3832 while( 1 ){
3833 ptr = my_fgets( buf , LBUF , fts ) ; /* read */
3834 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
3835
3836 fvec = decode_linebuf( buf ) ;
3837 if( fvec == NULL ) break ;
3838 if( fvec->nar == 0 ){ KILL_floatvec(fvec); break; }
3839
3840 for( ii=0 ; ii < fvec->nar && ii < ncol ; ii++ )
3841 tsar[nrow*ncol+ii] = fvec->ar[ii] ;
3842 for( ; ii < ncol ; ii++ )
3843 tsar[nrow*ncol+ii] = filler ; /* fill for incomplete lines */
3844 KILL_floatvec(fvec) ;
3845 nrow++ ; /* got one more complete row! */
3846 }
3847 fclose( fts ) ; /* finished with this file! */
3848 (void) my_fgets( NULL , 0 , NULL ) ; /* reset */
3849
3850 mri_add_name( fname , outim ) ;
3851 FRB(buf) ; lbfill = 0.0f ; RETURN(outim) ;
3852 }
3853
3854 /*---------------------------------------------------------------------------*/
3855 /*! Decode pairs of numbers separated by a single non-space character */
3856
decode_complex(char * str,float filler)3857 static INLINE complex decode_complex( char *str , float filler )
3858 {
3859 complex pp ; char ss ; float aa , bb ;
3860
3861 pp.r = pp.i = filler ;
3862 if( str == NULL ) return pp ;
3863 aa = bb = filler ;
3864 sscanf( str , "%f%c%f" , &aa , &ss , &bb ) ;
3865 pp.r = aa ; pp.i = bb ; return pp ;
3866 }
3867
3868 /*---------------------------------------------------------------------------*/
3869 /*! Ragged read pairs of values into a complex image. [08 Mar 2007] */
3870
mri_read_ascii_ragged_complex(char * fname,float filler)3871 MRI_IMAGE * mri_read_ascii_ragged_complex( char *fname , float filler )
3872 {
3873 MRI_IMAGE *outim ;
3874 complex *cxar , cval ;
3875 int ii,jj , ncol,nrow ;
3876 FILE *fts ;
3877 char *buf , *ptr ;
3878 NI_str_array *sar ; int nsar ;
3879
3880 ENTRY("mri_read_ascii_ragged_complex") ;
3881
3882 if( fname == NULL || *fname == '\0' ) RETURN(NULL) ;
3883
3884 fts = fopen(fname,"r");
3885 if( fts == NULL ){
3886 ERROR_message("mri_read_ascii_ragged_complex: couldn't open file %s" , fname ) ; RETURN(NULL) ;
3887 }
3888
3889 buf = (char *)malloc(LBUF) ;
3890
3891 /** step 1: read in ALL lines, see how many numbers are in each,
3892 in order to get the maximum row length and # of rows **/
3893
3894 (void) my_fgets( NULL , 0 , NULL ) ; /* reset */
3895 ncol = nrow = 0 ;
3896 while(1){
3897 ptr = my_fgets( buf , LBUF , fts ) ; /* read line */
3898 if( ptr==NULL || *ptr=='\0' ) break ; /* fails? end of data */
3899 sar = NI_decode_string_list( buf , "~" ) ; /* break into pieces */
3900 if( sar != NULL ){
3901 nsar = sar->num ; /* number of pieces */
3902 if( nsar > 0 ){ nrow++; ncol = MAX(ncol,nsar); }
3903 NI_delete_str_array(sar) ; /* recycle this */
3904 }
3905 }
3906 if( nrow == 0 || ncol == 0 ){ fclose(fts); free(buf); RETURN(NULL); }
3907
3908 /** At this point, ncol is the number of pairs to be read from each line **/
3909
3910 rewind(fts) ; /* start over at top of file */
3911
3912 outim = mri_new( ncol , nrow , MRI_complex ) ;
3913 cxar = MRI_COMPLEX_PTR(outim) ;
3914
3915 /** read lines, convert to floats, store **/
3916
3917 nrow = 0 ; cval.r = cval.i = filler ;
3918 while( 1 ){
3919 ptr = my_fgets( buf , LBUF , fts ) ; /* read line */
3920 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
3921 sar = NI_decode_string_list( buf , "~" ) ; /* break up */
3922 if( sar != NULL ){
3923 nsar = sar->num ; /* number of pieces */
3924 for( ii=0 ; ii < nsar ; ii++ ) /* decode each piece */
3925 cxar[nrow*ncol+ii] = decode_complex( sar->str[ii] , filler ) ;
3926 for( ; ii < ncol ; ii++ )
3927 cxar[nrow*ncol+ii] = cval ; /* fill row with junk */
3928 NI_delete_str_array(sar) ; /* done with this */
3929 }
3930 nrow++ ; /* added one complete row */
3931 }
3932
3933 free(buf); fclose( fts ); (void) my_fgets(NULL,0,NULL); /* cleanup */
3934
3935 mri_add_name(fname,outim) ; RETURN(outim) ;
3936 }
3937
3938 /*---------------------------------------------------------------------------*/
3939 /*! Decode vectors of numbers separated by a single non-space character;
3940 return value is number of values actually decoded.
3941 vec==NULL is OK for testing (then no values are assigned to it, duh).
3942 *//*-------------------------------------------------------------------------*/
3943
decode_fvect(char * str,float filler,int vdim,float * vec)3944 static int decode_fvect( char *str, float filler, int vdim, float *vec )
3945 {
3946 int ii,nn,mm ; float aa ;
3947
3948 if( vec != NULL ) for( ii=0 ; ii < vdim ; ii++ ) vec[ii] = filler ;
3949 if( str == NULL || *str == '\0' ) return 0 ;
3950
3951 if( *str == '*' ){ /* 23 Dec 2008 */
3952 ii = 1 ;
3953 if( str[1] == '*' && isdigit(str[2]) ) /* 16 Dec 2010 */
3954 ii = (int)strtod(str+2,NULL) ;
3955 return ii ;
3956 }
3957
3958 for( ii=0 ; ii < vdim ; ii++ ){
3959 nn = 0 ; mm = sscanf( str , "%f%n" , &aa , &nn ) ;
3960 if( mm == 0 ) return (ii) ;
3961 if( vec != NULL ) vec[ii] = aa ;
3962 str += nn ; if( *str == '\0' ) return (ii+1) ;
3963 str++ ; if( *str == '\0' ) return (ii+1) ;
3964 }
3965 return vdim ;
3966 }
3967
3968 /*---------------------------------------------------------------------------*/
3969 /*! Ragged read tuples of values into a fvect image.
3970 vdim = length of vectors to be read;
3971 can be zero, in which case will be determined from the data.
3972 *//*-------------------------------------------------------------------------*/
3973
mri_read_ascii_ragged_fvect(char * fname,float filler,int vdim)3974 MRI_IMAGE * mri_read_ascii_ragged_fvect( char *fname, float filler, int vdim )
3975 {
3976 MRI_IMAGE *outim ; float *var ;
3977 int ii,jj , ncol,nrow ;
3978 FILE *fts ;
3979 char *buf , *ptr ;
3980 NI_str_array *sar ; int nsar , nvdim ;
3981
3982 ENTRY("mri_read_ascii_ragged_fvect") ;
3983
3984 if( fname == NULL || *fname == '\0' ) RETURN(NULL) ;
3985
3986 if( strncmp(fname,"1D:",3) == 0 ){ /* cheap hack for 3dDeconvolve -stim_times */
3987 outim = mri_read_ragged_fromstring( fname+3 , filler ) ;
3988 if( outim != NULL && outim->kind == MRI_float){
3989 outim->kind = MRI_fvect ; outim->vdim = 1 ;
3990 }
3991 RETURN(outim) ;
3992 }
3993
3994 fts = fopen(fname,"r");
3995 if( fts == NULL ){
3996 ERROR_message("mri_read_ascii_ragged_fvect: couldn't open file %s" , fname ) ; RETURN(NULL) ;
3997 }
3998
3999 buf = (char *)malloc(LBUF) ;
4000
4001 /** step 1: read in ALL lines, see how many numbers are in each,
4002 in order to get the maximum row length and # of rows **/
4003
4004 (void) my_fgets( NULL , 0 , NULL ) ; /* reset */
4005 ncol = nrow = 0 ; nvdim = 0 ;
4006 while(1){
4007 ptr = my_fgets( buf , LBUF , fts ) ; /* read line */
4008 if( ptr==NULL || *ptr=='\0' ) break ; /* fails? end of data */
4009 sar = NI_decode_string_list( buf , "~" ) ; /* break into pieces */
4010 if( sar != NULL ){
4011 nsar = sar->num ; /* number of pieces */
4012 if( nsar > 0 ){ nrow++; ncol = MAX(ncol,nsar); }
4013 if( nsar > 0 && vdim == 0 ){ /* number of components */
4014 for( jj=0 ; jj < nsar ; jj++ ){
4015 ii = decode_fvect( sar->str[jj] , filler , 9999 , NULL ) ;
4016 nvdim = MAX(nvdim,ii) ;
4017 }
4018 }
4019 NI_delete_str_array(sar) ; /* recycle this */
4020 }
4021 }
4022 if( vdim == 0 ) vdim = nvdim ; /* set vdim from data */
4023 if( nrow == 0 || ncol == 0 || vdim == 0 ){
4024 fclose(fts); free(buf); RETURN(NULL);
4025 }
4026
4027 /** At this point, ncol is the number of vectors to be read from each line **/
4028
4029 rewind(fts) ; /* start over at top of file */
4030
4031 outim = mri_new_fvectim( ncol , nrow , 1 , vdim ) ;
4032 var = (float *)outim->im ;
4033 for( ii=0 ; ii < ncol*nrow*vdim ; ii++ ) var[ii] = filler ;
4034
4035 /** read lines, convert to floats, store **/
4036
4037 nrow = 0 ;
4038 while( 1 ){
4039 ptr = my_fgets( buf , LBUF , fts ) ; /* read line */
4040 if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
4041 sar = NI_decode_string_list( buf , "~" ) ; /* break up */
4042 if( sar != NULL ){
4043 nsar = sar->num ; /* number of pieces */
4044 for( ii=0 ; ii < nsar ; ii++ ) /* decode each piece */
4045 (void)decode_fvect( sar->str[ii] , filler , vdim ,
4046 var + (nrow*ncol+ii)*vdim ) ;
4047 NI_delete_str_array(sar) ; /* done with this */
4048 }
4049 nrow++ ; /* added one complete row */
4050 }
4051
4052 free(buf); fclose( fts ); (void) my_fgets(NULL,0,NULL); /* cleanup */
4053
4054 mri_add_name(fname,outim) ; RETURN(outim) ;
4055 }
4056
4057 /*---------------------------------------------------------------------------
4058 Read in an ASCII file to a float array.
4059 -----------------------------------------------------------------------------*/
4060
read_ascii_floats(char * fname,int * nff,float ** ff)4061 static void read_ascii_floats( char * fname, int * nff , float ** ff )
4062 {
4063 int ii,jj,val , used_tsar , alloc_tsar ;
4064 float *tsar ;
4065 float ftemp ;
4066 FILE *fts ;
4067 char *buf ; /* 08 Jul 2004: malloc this now, instead of auto */
4068 char *ptr ;
4069 int bpos , blen , nrow ;
4070
4071 /* check inputs */
4072
4073 if( nff == NULL || ff == NULL ) return ;
4074 if( fname == NULL || fname[0] == '\0' ){ *nff=0 ; *ff=NULL ; return ; }
4075
4076 fts = fopen( fname , "r" ) ;
4077 if( fts == NULL ){ *nff=0 ; *ff=NULL ; return ; }
4078 if( fts == NULL ){
4079 *nff = 0 ; *ff = NULL ;
4080 ERROR_message("mri_read_ascii_floats: couldn't open file %s" , fname ) ; return ;
4081 }
4082
4083 /* make some space */
4084
4085 used_tsar = 0 ;
4086 alloc_tsar = INC_TSARSIZE ;
4087 tsar = (float *) malloc( sizeof(float) * alloc_tsar ) ;
4088 if( tsar == NULL ){
4089 fprintf(stderr,"\n*** malloc fails: read_ascii_floats ***\n"); EXIT(1);
4090 }
4091
4092 /** read lines, convert to floats, store **/
4093
4094 nrow = 0 ;
4095 buf = (char *)malloc(LBUF) ;
4096 while( 1 ){
4097 ptr = afni_fgets( buf , LBUF , fts ) ; /* read */
4098 if( ptr == NULL ) break ; /* failure --> end of data */
4099 blen = strlen(buf) ;
4100 if( blen <= 0 ) break ; /* nothing --> end of data */
4101
4102 for( ii=0 ; ii < blen && isspace(buf[ii]) ; ii++ ) ; /* skip blanks */
4103
4104 if( ii == blen ) continue ; /* skip all blank line */
4105 if( buf[ii] == '#' ) continue ; /* skip a comment line */
4106 if( buf[ii] == '!' ) continue ;
4107
4108 /* convert commas to blanks */
4109
4110 for( jj=ii ; jj < blen ; jj++ ) if( buf[jj] == ',' ) buf[jj] = ' ' ;
4111
4112 for( bpos=ii ; bpos < blen ; ){
4113 val = sscanf( buf+bpos , "%f%n" , &ftemp , &jj ) ; /* read from string */
4114 if( val < 1 ) break ; /* bad read? */
4115 bpos += jj ; /* start of next read */
4116
4117 if( used_tsar == alloc_tsar ){
4118 alloc_tsar += INC_TSARSIZE ;
4119 tsar = (float *)realloc( tsar,sizeof(float)*alloc_tsar );
4120 if( tsar == NULL ){
4121 fprintf(stderr,"\n*** realloc fails: read_ascii_floats ***\n"); EXIT(1);
4122 }
4123 }
4124
4125 tsar[used_tsar++] = ftemp ; /* store input */
4126 }
4127
4128 nrow++ ; /* got one more complete row! */
4129 }
4130 fclose( fts ) ; /* finished with this file! */
4131 free( buf ) ;
4132
4133 if( used_tsar <= 1 ){ free(tsar); *nff=0; *ff=NULL; return; }
4134
4135 tsar = (float *) realloc( tsar , sizeof(float) * used_tsar ) ;
4136 if( tsar == NULL ){
4137 fprintf(stderr,"\n*** final realloc fails: read_ascii_floats ***\n"); EXIT(1);
4138 }
4139
4140 *nff = used_tsar; *ff = tsar; return;
4141 }
4142
4143 /*--------------------------------------------------------------
4144 Read a pile of images from one ASCII file.
4145 Adapted from mri_read_3D - Feb 2000 - RWCox.
4146 [N.B.: if this routine is altered, don't forget mri_imcount!]
4147 ----------------------------------------------------------------*/
4148
mri_read_3A(char * tname)4149 MRI_IMARR * mri_read_3A( char *tname )
4150 {
4151 int nx , ny , nz , ii , nxyz,nxy , nff ;
4152 int ngood , length , kim , datum_type ;
4153 char fname[256]="\0" , buf[512] ;
4154 MRI_IMARR *newar ;
4155 MRI_IMAGE *newim=NULL , * flim ;
4156 float * ff ;
4157
4158 ENTRY("mri_read_3A") ;
4159
4160 /*** get info from 3A tname ***/
4161
4162 if( tname == NULL || strlen(tname) < 10 ) RETURN(NULL) ;
4163
4164 switch( tname[2] ){ /* allow for 3As:, 3Ab:, 3Af: */
4165
4166 default: ngood = 0 ; break ;
4167
4168 case 's':
4169 ngood = sscanf( tname, "3As:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
4170 datum_type = MRI_short ;
4171 break ;
4172
4173 case 'b':
4174 ngood = sscanf( tname, "3Ab:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
4175 datum_type = MRI_byte ;
4176 break ;
4177
4178 case 'f':
4179 ngood = sscanf( tname, "3Af:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
4180 datum_type = MRI_float ;
4181 break ;
4182 }
4183
4184 if( ngood < 4 || nx <= 0 || ny <= 0 || nz <= 0 || strlen(fname) <= 0 ) RETURN(NULL) ;
4185
4186 /* read the input file */
4187
4188 read_ascii_floats( fname , &nff , &ff ) ;
4189
4190 if( nff <= 0 || ff == NULL ) RETURN(NULL) ;
4191
4192 nxy = nx*ny ; nxyz = nxy*nz ;
4193
4194 if( nff < nxyz ){
4195 fprintf(stderr,
4196 "\n** WARNING: %s is too short - padding with %d zeros\n",
4197 tname,nxyz-nff) ;
4198 ff = (float *) realloc( ff , sizeof(float) * nxyz ) ;
4199 for( ii=nff ; ii < nxyz ; ii++ ) ff[ii] = 0.0 ;
4200 nff = nxyz ;
4201 } else if( nff > nxyz ){
4202 fprintf(stderr,
4203 "\n** WARNING: %s is too long - truncating off last %d values\n",
4204 tname,nff-nxyz) ;
4205 }
4206
4207 /* put the input data into MRI_IMAGEs */
4208
4209 INIT_IMARR(newar) ;
4210
4211 for( kim=0 ; kim < nz ; kim++ ){
4212 flim = mri_new( nx,ny , MRI_float ) ;
4213 memcpy( MRI_FLOAT_PTR(flim) , ff+nxy*kim , sizeof(float)*nxy ) ;
4214 switch( datum_type ){
4215 case MRI_float: newim = flim ; break ;
4216 case MRI_short: newim = mri_to_short(1.0,flim) ; mri_free(flim) ; break ;
4217 case MRI_byte: newim = mri_to_byte_scl(1.0,0.0,flim) ; mri_free(flim) ; break ;
4218 }
4219
4220 if( nz == 1 ) mri_add_name( fname , newim ) ;
4221 else {
4222 sprintf( buf , "%s#%d" , fname,kim ) ;
4223 mri_add_name( buf , newim ) ;
4224 }
4225
4226 ADDTO_IMARR(newar,newim) ;
4227 }
4228
4229 free(ff) ; RETURN(newar) ;
4230 }
4231
4232 /*---------------------------------------------------------------------------
4233 Stuff to read an ANALYZE 7.5 .img file, given the .hdr filename
4234 -- 05 Feb 2001 - RWCox
4235 -- 27 Nov 2001 - modified to use funused1 as a scale factor
4236 -----------------------------------------------------------------------------*/
4237
4238 #include "mayo_analyze.h"
4239
4240 /*---------------------------------------------------------------*/
4241 /*! Byte swap ANALYZE file header in various places */
4242
swap_analyze_hdr(struct dsr * pntr)4243 static void swap_analyze_hdr( struct dsr *pntr )
4244 {
4245 ENTRY("swap_analyze_hdr") ;
4246 swap_4(&pntr->hk.sizeof_hdr) ;
4247 swap_4(&pntr->hk.extents) ;
4248 swap_2(&pntr->hk.session_error) ;
4249 swap_2(&pntr->dime.dim[0]) ;
4250 swap_2(&pntr->dime.dim[1]) ;
4251 swap_2(&pntr->dime.dim[2]) ;
4252 swap_2(&pntr->dime.dim[3]) ;
4253 swap_2(&pntr->dime.dim[4]) ;
4254 swap_2(&pntr->dime.dim[5]) ;
4255 swap_2(&pntr->dime.dim[6]) ;
4256 swap_2(&pntr->dime.dim[7]) ;
4257 #if 0
4258 swap_2(&pntr->dime.unused1) ;
4259 #endif
4260 swap_2(&pntr->dime.datatype) ;
4261 swap_2(&pntr->dime.bitpix) ;
4262 swap_4(&pntr->dime.pixdim[0]) ;
4263 swap_4(&pntr->dime.pixdim[1]) ;
4264 swap_4(&pntr->dime.pixdim[2]) ;
4265 swap_4(&pntr->dime.pixdim[3]) ;
4266 swap_4(&pntr->dime.pixdim[4]) ;
4267 swap_4(&pntr->dime.pixdim[5]) ;
4268 swap_4(&pntr->dime.pixdim[6]) ;
4269 swap_4(&pntr->dime.pixdim[7]) ;
4270 swap_4(&pntr->dime.vox_offset) ;
4271 swap_4(&pntr->dime.funused1) ;
4272 swap_4(&pntr->dime.funused2) ;
4273 swap_4(&pntr->dime.cal_max) ;
4274 swap_4(&pntr->dime.cal_min) ;
4275 swap_4(&pntr->dime.compressed) ;
4276 swap_4(&pntr->dime.verified) ;
4277 swap_2(&pntr->dime.dim_un0) ;
4278 swap_4(&pntr->dime.glmax) ;
4279 swap_4(&pntr->dime.glmin) ;
4280 EXRETURN ;
4281 }
4282
4283 /*---------------------------------------------------------------*/
4284 /*! Count how many 2D slices are in an ANALYZE file.
4285
4286 \param hname = the "hdr" file of the hdr/img file pair.
4287 */
4288
mri_imcount_analyze75(char * hname)4289 static int mri_imcount_analyze75( char * hname )
4290 {
4291 FILE * fp ;
4292 struct dsr hdr ; /* ANALYZE .hdr format */
4293 int doswap , nz ;
4294
4295 ENTRY("mri_imcount_analyze75") ;
4296
4297 fp = fopen( hname , "rb" ) ;
4298 if( fp == NULL ){
4299 ERROR_message("mri_imcount_analyze75: couldn't open file %s" , hname ) ; RETURN(0) ;
4300 }
4301 hdr.dime.dim[0] = 0 ;
4302 fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
4303 fclose(fp) ;
4304 if( hdr.dime.dim[0] == 0 ) RETURN(0) ;
4305 doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
4306 if( doswap ) swap_analyze_hdr( &hdr ) ;
4307
4308 switch( hdr.dime.dim[0] ){
4309 case 2: nz = 1 ; break ;
4310 case 3: nz = hdr.dime.dim[3] ; break ;
4311
4312 default:
4313 case 4: nz = hdr.dime.dim[3] * hdr.dime.dim[4] ; break ;
4314 }
4315 if( nz < 1 ) nz = 1 ;
4316
4317 RETURN(nz) ;
4318 }
4319
4320 /*---------------------------------------------------------------*/
4321 /*! Read an ANALYZE file into an ARRAY of 2D images.
4322
4323 \param hname = the "hdr" file for the hdr/img pair
4324 */
4325
mri_read_analyze75(char * hname)4326 MRI_IMARR * mri_read_analyze75( char * hname )
4327 {
4328 FILE * fp ;
4329 char iname[1024] , buf[1024] ;
4330 int ii , jj , doswap ;
4331 struct dsr hdr ; /* ANALYZE .hdr format */
4332 int ngood , length , kim , koff , datum_type , datum_len , swap ;
4333 int nx,ny,nz , hglobal=0 , himage=0 ;
4334 float dx,dy,dz ;
4335 MRI_IMARR * newar ;
4336 MRI_IMAGE * newim ;
4337 void * imar ;
4338 float fac=0.0 ; /* 27 Nov 2001 */
4339 int floatize ; /* 28 Nov 2001 */
4340 int spmorg=0 ; /* 28 Nov 2001 */
4341
4342 ENTRY("mri_read_analyze75") ;
4343
4344 /* check & prepare filenames */
4345
4346 if( hname == NULL ) RETURN(NULL) ;
4347 jj = strlen(hname) ;
4348 if( jj < 5 ) RETURN(NULL) ;
4349 if( strcmp(hname+jj-3,"hdr") != 0 ) RETURN(NULL) ;
4350 strcpy(iname,hname) ; strcpy(iname+jj-3,"img") ;
4351
4352 /* read header file into struct */
4353
4354 fp = fopen( hname , "rb" ) ;
4355 if( fp == NULL ){
4356 ERROR_message("mri_read_analyze75: couldn't open file %s" , hname ) ; RETURN(NULL) ;
4357 }
4358 hdr.dime.dim[0] = 0 ;
4359 fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
4360 fclose(fp) ;
4361 if( hdr.dime.dim[0] == 0 ) RETURN(NULL) ;
4362
4363 /* check for swap-age */
4364
4365 doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
4366 if( doswap ) swap_analyze_hdr( &hdr ) ;
4367
4368 /* 28 Nov 2001: attempt to decode originator a la SPM */
4369
4370 { short xyzuv[5] , xx,yy,zz ;
4371 memcpy( xyzuv , hdr.hist.originator , 10 ) ;
4372 if( xyzuv[3] == 0 && xyzuv[4] == 0 ){
4373 xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
4374 if( doswap ){ swap_2(&xx); swap_2(&yy); swap_2(&zz); }
4375 if( xx > 0 && xx < hdr.dime.dim[1] &&
4376 yy > 0 && yy < hdr.dime.dim[2] &&
4377 zz > 0 && zz < hdr.dime.dim[3] ) spmorg = 1 ;
4378 }
4379 }
4380 if( spmorg ) strcpy( MRILIB_orients , "LRPAIS" ) ;
4381
4382 /* 27 Nov 2001: get a scale factor for images */
4383
4384 if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
4385 fac = hdr.dime.funused1 ;
4386 (void) thd_floatscan( 1 , &fac ) ;
4387 if( fac < 0.0f || fac == 1.0f ) fac = 0.0f ;
4388 else if( fac < 1.e-6f || fac > 1.e+6f ){
4389 WARNING_message("ANALYZE file %s has scale factor (funused1) set to %g",hname,fac) ;
4390 ININFO_message ("to ignore this scale factor, setenv AFNI_ANALYZE_SCALE NO") ;
4391 }
4392 }
4393
4394 floatize = (fac != 0.0 || AFNI_yesenv("AFNI_ANALYZE_FLOATIZE")) ; /* 28 Nov 2001 */
4395
4396 /* get data type into mrilib MRI_* form */
4397
4398 switch( hdr.dime.datatype ){
4399 default:
4400 fprintf(stderr,"*** %s: Unknown ANALYZE datatype=%d (%s)\n",
4401 hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
4402 RETURN(NULL) ;
4403
4404 case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte ; break;
4405 case ANDT_SIGNED_SHORT: datum_type = MRI_short ; break;
4406 case ANDT_SIGNED_INT: datum_type = MRI_int ; break;
4407 case ANDT_FLOAT: datum_type = MRI_float ; floatize = 0; break;
4408 case ANDT_COMPLEX: datum_type = MRI_complex; floatize = 0; break;
4409 case ANDT_RGB: datum_type = MRI_rgb ; floatize = 0; break;
4410 case ANDT_DOUBLE: datum_type = MRI_double ; floatize = 0; break;
4411 }
4412
4413 datum_len = mri_datum_size(datum_type) ;
4414
4415 /* compute dimensions of images, and number of images */
4416
4417 nx = hdr.dime.dim[1] ;
4418 ny = hdr.dime.dim[2] ;
4419 if( nx < 2 || ny < 2 ) RETURN(NULL) ;
4420
4421 switch( hdr.dime.dim[0] ){
4422 case 2: nz = 1 ; break ;
4423 case 3: nz = hdr.dime.dim[3] ; break ;
4424
4425 default:
4426 case 4: nz = hdr.dime.dim[3] * hdr.dime.dim[4] ; break ;
4427 }
4428 if( nz < 1 ) nz = 1 ;
4429
4430 dx = hdr.dime.pixdim[1] ;
4431 dy = hdr.dime.pixdim[2] ;
4432 dz = hdr.dime.pixdim[3] ;
4433
4434 /** fprintf(stderr,"mri_read_analyze75: nx=%d ny=%d nz=%d\n",nx,ny,nz) ; **/
4435 /** fprintf(stderr,"mri_read_analyze75: dx=%g dy=%g dz=%g\n",dx,dy,dz) ; **/
4436
4437 /* open .img file and read images from it */
4438
4439 length = THD_filesize(iname) ;
4440 if( length <= 0 ){
4441 fprintf(stderr,"*** Can't find ANALYZE file %s\n",iname) ;
4442 RETURN(NULL) ;
4443 }
4444
4445 fp = fopen( iname , "rb" ) ;
4446 if( fp == NULL ){
4447 ERROR_message("mri_read_analyze75: couldn't open file %s" , iname ) ; RETURN(NULL) ;
4448 }
4449
4450 ngood = datum_len*nx*ny*nz ;
4451 if( length < ngood ){
4452 fprintf( stderr,
4453 "*** ANALYZE file %s is %d bytes long but must be at least %d bytes long\n"
4454 "*** for nx=%d ny=%d nz=%d and voxel=%d bytes\n",
4455 iname,length,ngood,nx,ny,nz,datum_len ) ;
4456 fclose(fp) ; RETURN(NULL) ;
4457 }
4458
4459 /*** read images from the file ***/
4460
4461 INIT_IMARR(newar) ;
4462
4463 for( kim=0 ; kim < nz ; kim++ ){
4464 koff = hglobal + (kim+1)*himage + datum_len*nx*ny*kim ;
4465 /** fprintf(stderr,"mri_read_analyze75: kim=%d koff=%d\n",kim,koff) ; **/
4466 fseek( fp , koff , SEEK_SET ) ;
4467
4468 newim = mri_new( nx , ny , datum_type ) ;
4469 imar = mri_data_pointer( newim ) ;
4470 length = fread( imar , datum_len , nx * ny , fp ) ;
4471
4472 if( doswap ){
4473 switch( datum_len ){
4474 default: break ;
4475 case 2: swap_twobytes ( nx*ny , imar ) ; break ; /* short */
4476 case 4: swap_fourbytes( nx*ny , imar ) ; break ; /* int, float */
4477 case 8: swap_fourbytes( 2*nx*ny , imar ) ; break ; /* complex */
4478 }
4479 newim->was_swapped = 1 ; /* 07 Mar 2002 */
4480 }
4481
4482 /* 28 Nov 2001: convert to floats? */
4483
4484 if( floatize ){
4485 MRI_IMAGE *qim = mri_to_float(newim) ;
4486 mri_free(newim) ; newim = qim ;
4487 }
4488
4489 if( nz == 1 ) mri_add_name( iname , newim ) ;
4490 else {
4491 sprintf( buf , "%s#%d" , iname,kim ) ;
4492 mri_add_name( buf , newim ) ;
4493 }
4494
4495 newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dw = 1.0 ;
4496 ADDTO_IMARR(newar,newim) ;
4497
4498 /* 27 Nov 2001: scale image? */
4499
4500 if( fac != 0.0 ) mri_scale_inplace( fac , newim ) ;
4501 }
4502
4503 fclose(fp) ; RETURN(newar) ;
4504 }
4505
4506 /*-----------------------------------------------------------------*/
4507 /*! Read an ANALYZE file into an ARRAY of 3D images [26 Aug 2002].
4508
4509 \param hname = the "hdr" file for the hdr/img pair
4510 */
4511
mri_read3D_analyze75(char * hname)4512 MRI_IMARR * mri_read3D_analyze75( char * hname )
4513 {
4514 FILE * fp ;
4515 char iname[1024] , buf[1024] ;
4516 int ii , jj , doswap ;
4517 struct dsr hdr ; /* ANALYZE .hdr format */
4518 int ngood , length , kim , koff , datum_type , datum_len , swap ;
4519 int nx,ny,nz , hglobal=0 , himage=0 ;
4520 float dx,dy,dz ;
4521 MRI_IMARR * newar ;
4522 MRI_IMAGE * newim ;
4523 void * imar ;
4524 float fac=0.0 ; /* 27 Nov 2001 */
4525 int floatize ; /* 28 Nov 2001 */
4526 int spmorg=0 ; /* 28 Nov 2001 */
4527
4528 int nt , nxyz ; /* 26 Aug 2002 */
4529 float dt ;
4530
4531 ENTRY("mri_read3D_analyze75") ;
4532
4533 /* check & prepare filenames */
4534
4535 if( hname == NULL ) RETURN(NULL) ;
4536 jj = strlen(hname) ;
4537 if( jj < 5 ) RETURN(NULL) ;
4538 if( strcmp(hname+jj-3,"hdr") != 0 ) RETURN(NULL) ;
4539 strcpy(iname,hname) ; strcpy(iname+jj-3,"img") ;
4540
4541 /* read header file into struct */
4542
4543 fp = fopen( hname , "rb" ) ;
4544 if( fp == NULL ){
4545 ERROR_message("mri_read3D_analyze75: couldn't open file %s" , hname ) ; RETURN(NULL) ;
4546 }
4547 hdr.dime.dim[0] = 0 ;
4548 fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
4549 fclose(fp) ;
4550 if( hdr.dime.dim[0] == 0 ) RETURN(NULL) ;
4551
4552 /* check for swap-age */
4553
4554 doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
4555 if( doswap ) swap_analyze_hdr( &hdr ) ;
4556
4557 /* 28 Nov 2001: attempt to decode originator a la SPM */
4558
4559 { short xyzuv[5] , xx,yy,zz ;
4560 memcpy( xyzuv , hdr.hist.originator , 10 ) ;
4561 if( xyzuv[3] == 0 && xyzuv[4] == 0 ){
4562 xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
4563 if( doswap ){ swap_2(&xx); swap_2(&yy); swap_2(&zz); }
4564 if( xx > 0 && xx < hdr.dime.dim[1] &&
4565 yy > 0 && yy < hdr.dime.dim[2] &&
4566 zz > 0 && zz < hdr.dime.dim[3] ) spmorg = 1 ;
4567 }
4568 }
4569 if( spmorg ) strcpy( MRILIB_orients , "LRPAIS" ) ;
4570
4571 /* 27 Nov 2001: get a scale factor for images */
4572
4573 if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
4574 fac = hdr.dime.funused1 ;
4575 (void) thd_floatscan( 1 , &fac ) ;
4576 if( fac < 0.0 || fac == 1.0 ) fac = 0.0 ;
4577 }
4578
4579 floatize = (fac != 0.0 || AFNI_yesenv("AFNI_ANALYZE_FLOATIZE")) ; /* 28 Nov 2001 */
4580
4581 /* get data type into mrilib MRI_* form */
4582
4583 switch( hdr.dime.datatype ){
4584 default:
4585 fprintf(stderr,"*** %s: Unknown ANALYZE datatype=%d (%s)\n",
4586 hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
4587 RETURN(NULL) ;
4588
4589 case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte ; break;
4590 case ANDT_SIGNED_SHORT: datum_type = MRI_short ; break;
4591 case ANDT_SIGNED_INT: datum_type = MRI_int ; break;
4592 case ANDT_FLOAT: datum_type = MRI_float ; floatize = 0; break;
4593 case ANDT_COMPLEX: datum_type = MRI_complex; floatize = 0; break;
4594 case ANDT_RGB: datum_type = MRI_rgb ; floatize = 0; break;
4595 }
4596
4597 datum_len = mri_datum_size(datum_type) ;
4598
4599 /* compute dimensions of images, and number of images */
4600
4601 nx = hdr.dime.dim[1] ;
4602 ny = hdr.dime.dim[2] ;
4603 if( nx < 2 || ny < 2 ) RETURN(NULL) ;
4604
4605 switch( hdr.dime.dim[0] ){
4606 case 2: nz = 1 ; nt = 1 ; ; break ;
4607 case 3: nz = hdr.dime.dim[3] ; nt = 1 ; ; break ;
4608
4609 default:
4610 case 4: nz = hdr.dime.dim[3] ; nt = hdr.dime.dim[4] ; break ;
4611 }
4612 if( nz < 1 ) nz = 1 ;
4613 if( nt < 1 ) nt = 1 ;
4614
4615 dx = hdr.dime.pixdim[1] ;
4616 dy = hdr.dime.pixdim[2] ;
4617 dz = hdr.dime.pixdim[3] ;
4618 dt = hdr.dime.pixdim[4] ; if( dt <= 0.0 ) dt = 1.0 ;
4619
4620 /* open .img file and read images from it */
4621
4622 length = THD_filesize(iname) ;
4623 if( length <= 0 ){
4624 fprintf(stderr,"*** Can't find ANALYZE file %s\n",iname) ;
4625 RETURN(NULL) ;
4626 }
4627
4628 fp = fopen( iname , "rb" ) ;
4629 if( fp == NULL ){
4630 ERROR_message("mri_read3D_analyze75: couldn't open file %s" , iname ) ; RETURN(NULL) ;
4631 }
4632
4633 ngood = datum_len*nx*ny*nz*nt ;
4634 if( length < ngood ){
4635 fprintf( stderr,
4636 "*** ANALYZE file %s is %d bytes long but must be at least %d bytes long\n"
4637 "*** for nx=%d ny=%d nz=%d nt=%d and voxel=%d bytes\n",
4638 iname,length,ngood,nx,ny,nz,nt,datum_len ) ;
4639 fclose(fp) ; RETURN(NULL) ;
4640 }
4641
4642 /*** read images from the file ***/
4643
4644 INIT_IMARR(newar) ;
4645
4646 nxyz = nx*ny*nz ;
4647 for( kim=0 ; kim < nt ; kim++ ){
4648 koff = hglobal + (kim+1)*himage + datum_len*nxyz*kim ;
4649 fseek( fp , koff , SEEK_SET ) ;
4650
4651 newim = mri_new_vol( nx,ny,nz , datum_type ) ;
4652 imar = mri_data_pointer( newim ) ;
4653 length = fread( imar , datum_len , nxyz , fp ) ;
4654
4655 if( doswap ){
4656 switch( datum_len ){
4657 default: break ;
4658 case 2: swap_twobytes ( nxyz , imar ) ; break ; /* short */
4659 case 4: swap_fourbytes( nxyz , imar ) ; break ; /* int, float */
4660 case 8: swap_fourbytes( 2*nxyz , imar ) ; break ; /* complex */
4661 }
4662 newim->was_swapped = 1 ; /* 07 Mar 2002 */
4663 }
4664
4665 /* 28 Nov 2001: convert to floats? */
4666
4667 if( floatize ){
4668 MRI_IMAGE *qim = mri_to_float(newim) ;
4669 mri_free(newim) ; newim = qim ;
4670 }
4671
4672 if( nt == 1 ) mri_add_name( iname , newim ) ;
4673 else {
4674 sprintf( buf , "%s#%d" , iname,kim ) ;
4675 mri_add_name( buf , newim ) ;
4676 }
4677
4678 newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dt = dt ; newim->dw = 1.0 ;
4679 ADDTO_IMARR(newar,newim) ;
4680
4681 /* 27 Nov 2001: scale image? */
4682
4683 if( fac != 0.0 ) mri_scale_inplace( fac , newim ) ;
4684 }
4685
4686 fclose(fp) ; RETURN(newar) ;
4687 }
4688
4689 /*---------------------------------------------------------------------------
4690 12 Mar 2001 - stuff to read a Siemens Vision .ima file
4691 -----------------------------------------------------------------------------*/
4692
4693 #include "siemens_vision.h"
4694
4695 /*! Count the number of 2D images in a Siemens Vision .ima file.
4696
4697 Unfortunately, this requires reading the image data and checking
4698 for all-zero images. This is because Siemens stores their data
4699 in a fixed size file, and so just fills out the empty space with
4700 blank images if need be.
4701 */
4702
mri_imcount_siemens(char * hname)4703 static int mri_imcount_siemens( char * hname )
4704 {
4705 struct Siemens_vision_header head ;
4706 FILE * fp ;
4707 int i,j,xx,yy , matrix , swap , imagesize,nxx,blank , slices ;
4708 struct stat file_stat ;
4709 short *imar ;
4710
4711 /*--- check file size ---*/
4712
4713 if( hname == NULL ) return 0 ;
4714
4715 i = stat( hname , &file_stat ) ;
4716 if( i < 0 ) return 0 ;
4717
4718 /*--- read header data ---*/
4719
4720 fp = fopen( hname , "rb" ) ;
4721 if( fp == NULL ){
4722 ERROR_message("mri_imcount_siemens: couldn't open file %s" , hname ) ; return 0 ;
4723 }
4724 fread( &head , sizeof(struct Siemens_vision_header) , 1 , fp ) ;
4725
4726 /*-- check some integer in header to determine if we need to byteswap --*/
4727
4728 swap = ( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ) ;
4729 if( swap ){
4730 swap_4( &(head.SiemensStudyDateMM) ) ;
4731 if( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ){
4732 swap = 0 ;
4733 }
4734 }
4735
4736 /*-- find image size from header --*/
4737
4738 if( swap ) swap_4( &(head.DisplayMatrixSize) ) ;
4739 imagesize = head.DisplayMatrixSize ;
4740
4741 /*-- determine number of sub-images in file --*/
4742
4743 #undef MATRIX_MAX
4744 #define MATRIX_MAX 9
4745
4746 i = 2*imagesize*imagesize ;
4747 for( matrix=1 ; matrix < MATRIX_MAX ; matrix++ )
4748 if( file_stat.st_size == i*matrix*matrix + SIEMENS_HEADERSIZE ) break ;
4749
4750 if( matrix == MATRIX_MAX ){
4751 fclose(fp) ; return 0 ; /* didn't recognize file format */
4752 }
4753 #undef MATRIX_MAX
4754
4755 /*-- read image data from file (but don't byteswap it) --*/
4756
4757 imar = (short *) calloc(sizeof(short),matrix*matrix*imagesize*imagesize) ;
4758 fseek( fp , SIEMENS_HEADERSIZE , SEEK_SET ) ;
4759 fread( imar , sizeof(short) , matrix*matrix*imagesize*imagesize , fp ) ;
4760 fclose(fp) ;
4761
4762 /*-- count slices - all zero (blank) slices at end are skipped --*/
4763
4764 slices = 0 ; nxx = matrix*imagesize ;
4765
4766 for( yy=0 ; yy < matrix ; yy++ ){ /* rows in array of sub-images */
4767 for( xx=0 ; xx < matrix ; xx++ ){ /* cols in array of sub-images */
4768 blank = 1 ;
4769 for( j=0 ; j < imagesize ; j++ ){ /* row in sub-image */
4770 for( i=0 ; i < imagesize ; i++ ){ /* col in sub-image */
4771 if( imar[i+xx*imagesize+(j+yy*imagesize)*nxx] ) blank = 0 ;
4772 }
4773 }
4774 if( !blank ) slices = 1 + xx + yy*matrix ;
4775 }
4776 }
4777
4778 free(imar) ; return slices ;
4779 }
4780
4781 /*---------------------------------------------------------------------------*/
4782 /*! Read an array of 2D images from Siemens Vision .ima file.
4783
4784 The images are stored in a 2D array, which requires untangling the
4785 data rows to put them into separate MRI_IMAGE structs.
4786 */
4787
mri_read_siemens(char * hname)4788 MRI_IMARR * mri_read_siemens( char * hname )
4789 {
4790 struct Siemens_vision_header head ;
4791 FILE * fp ;
4792 int i,j,xx,yy , matrix , swap , imagesize,nxx,blank , slices,nz ;
4793 struct stat file_stat ;
4794 short *imar ;
4795 MRI_IMARR * newar ;
4796 MRI_IMAGE * newim ;
4797 short * nar ;
4798 char buf[256] ;
4799 float dx,dy,dz ;
4800 char *eee ; int ileave=0 ; /* 25 Sep 2001 */
4801
4802 ENTRY("mri_read_siemens") ;
4803
4804 /*--- check file size ---*/
4805
4806 if( hname == NULL ) RETURN(NULL) ;
4807
4808 i = stat( hname , &file_stat ) ;
4809 if( i < 0 ) RETURN(NULL) ;
4810
4811 /*--- read header data ---*/
4812
4813 fp = fopen( hname , "rb" ) ;
4814 if( fp == NULL ){
4815 ERROR_message("mri_read_siemens: couldn't open file %s" , hname ) ; RETURN(NULL) ;
4816 }
4817 fread( &head , sizeof(struct Siemens_vision_header) , 1 , fp ) ;
4818
4819 /*-- check some integer in header to determine if we need to byteswap --*/
4820
4821 swap = ( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ) ;
4822 if( swap ){
4823 swap_4( &(head.SiemensStudyDateMM) ) ;
4824 if( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ){
4825 swap = 0 ;
4826 }
4827 }
4828
4829 /*-- find image size from header --*/
4830
4831 if( swap ) swap_4( &(head.DisplayMatrixSize) ) ;
4832 imagesize = head.DisplayMatrixSize ;
4833
4834 /*-- determine number of sub-images in file --*/
4835
4836 #undef MATRIX_MAX
4837 #define MATRIX_MAX 16
4838
4839 i = 2*imagesize*imagesize ;
4840 for( matrix=1 ; matrix < MATRIX_MAX ; matrix++ )
4841 if( file_stat.st_size == i*matrix*matrix + SIEMENS_HEADERSIZE ) break ;
4842
4843 if( matrix == MATRIX_MAX ){
4844 fclose(fp) ; RETURN(NULL) ; /* didn't recognize file format */
4845 }
4846 #undef MATRIX_MAX
4847
4848 /*-- read image data from file and byteswap it, if needed --*/
4849
4850 imar = (short *) calloc(sizeof(short),matrix*matrix*imagesize*imagesize) ;
4851 fseek( fp , SIEMENS_HEADERSIZE , SEEK_SET ) ;
4852 fread( imar , sizeof(short) , matrix*matrix*imagesize*imagesize , fp ) ;
4853 fclose(fp) ;
4854
4855 if( swap ) swap_twobytes( matrix*matrix*imagesize*imagesize , imar ) ;
4856
4857 /*-- count slices - all zero (blank) slices at end are skipped --*/
4858
4859 slices = 0 ; nxx = matrix*imagesize ;
4860
4861 for( yy=0 ; yy < matrix ; yy++ ){ /* rows in array of sub-images */
4862 for( xx=0 ; xx < matrix ; xx++ ){ /* cols in array of sub-images */
4863 blank = 1 ;
4864 for( j=0 ; j < imagesize ; j++ ){ /* row in sub-image */
4865 for( i=0 ; i < imagesize ; i++ ){ /* col in sub-image */
4866 if( imar[i+xx*imagesize+(j+yy*imagesize)*nxx] ) blank = 0 ;
4867 }
4868 }
4869 if( !blank ) slices = 1 + xx + yy*matrix ;
4870 }
4871 }
4872
4873 if( slices == 0 ){ free(imar) ; RETURN(NULL) ; } /* bad news */
4874
4875 /*-- get image dimensions, etc --*/
4876
4877 if( swap ){
4878 swap_8(&(head.FOVRow));
4879 swap_8(&(head.FOVColumn));
4880 swap_8(&(head.SliceThickness));
4881 }
4882 dx = head.FOVRow / imagesize ;
4883 dy = head.FOVColumn / imagesize ;
4884 dz = head.SliceThickness ;
4885
4886 /*-- save orientation and offset in global variables --*/
4887
4888 MRILIB_orients[0] = head.OrientationSet1Left[0] ;
4889 MRILIB_orients[1] = head.OrientationSet2Right[0];
4890 MRILIB_orients[2] = head.OrientationSet1Top[0] ;
4891 MRILIB_orients[3] = head.OrientationSet2Down[0] ;
4892 MRILIB_orients[4] = head.OrientationSet1Back[0] ;
4893 MRILIB_orients[5] = head.OrientationSet2Front[0];
4894 for (i=0; i<6; i++) {
4895 if (MRILIB_orients[i]=='H') MRILIB_orients[i]='S';
4896 if (MRILIB_orients[i]=='F') MRILIB_orients[i]='I';
4897 }
4898 MRILIB_orients[6] = '\0' ;
4899 MRILIB_zoff = fabs(strtod(head.TextSlicePosition,NULL)) ; use_MRILIB_zoff = 1 ;
4900
4901 /*-- create output --*/
4902
4903 INIT_IMARR(newar) ;
4904
4905 for( yy=0 ; yy < matrix ; yy++ ){ /* rows in array of sub-images */
4906 for( xx=0 ; xx < matrix ; xx++ ){ /* cols in array of sub-images */
4907
4908 newim = mri_new( imagesize , imagesize , MRI_short ) ;
4909 nar = MRI_SHORT_PTR( newim ) ;
4910
4911 if( swap ) newim->was_swapped = 1 ; /* 07 Mar 2002 */
4912
4913 for( j=0 ; j < imagesize ; j++ ) /* row in sub-image */
4914 memcpy( nar+j*imagesize ,
4915 imar+xx*imagesize+(j+yy*imagesize)*nxx , 2*imagesize ) ;
4916
4917 sprintf( buf , "%s#%d:%d" , hname,xx,yy ) ;
4918 mri_add_name( buf , newim ) ;
4919
4920 newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dw = 1.0 ;
4921 ADDTO_IMARR(newar,newim) ;
4922 if( IMARR_COUNT(newar) == slices ) goto Done ; /* Aauugghh!!! */
4923 }
4924 }
4925
4926 Done:
4927
4928 /*-- 25 Sep 2001: possibly interleave the images --*/
4929
4930 eee = getenv("AFNI_SIEMENS_INTERLEAVE") ;
4931 ileave = ( (eee != NULL) && (*eee=='Y' || *eee=='y') ) ;
4932 if( ileave && slices > 2 ){
4933 int mid = (slices-1)/2 ; /* midpoint */
4934 MRI_IMARR *qar ; /* new image array */
4935 INIT_IMARR(qar) ;
4936 for( i=0 ; i < slices ; i++ ){
4937 if( i%2 == 0 ) j = i/2 ; /* slice #i is in newar #j */
4938 else j = mid + (i+1)/2 ;
4939 ADDTO_IMARR(qar,IMARR_SUBIM(newar,j)) ; /* move image to new array */
4940 }
4941 FREE_IMARR(newar) ; newar = qar ;
4942 }
4943
4944 free(imar) ; RETURN(newar) ;
4945 }
4946
4947 /*---------------------------------------------------------------------------*/
4948 /*! Check for dicom magic number (string) in file
4949 Bytes 128-131 should be "DICM" in a Dicom Part 10 file
4950 */
4951
check_dicom_magic_num(char * fname)4952 static int check_dicom_magic_num( char *fname )
4953 {
4954 FILE *fp;
4955 char test_string[5] ;
4956
4957 fp = fopen( fname, "rb" ) ;
4958 if( fp == NULL ){
4959 ERROR_message("check_dicom_magic_num: couldn't open file %s" , fname ) ; return 0 ;
4960 }
4961 fseek( fp, 128 , SEEK_SET ) ;
4962 fread( test_string , 1 , 4 , fp ) ; test_string[4] = '\0' ;
4963 fclose( fp ) ;
4964 if( strcmp(test_string,"DICM") == 0 ) {
4965 return 1 ;
4966 } else {
4967 return 0 ;
4968 }
4969 }
4970
4971 /*---------------------------------------------------------------------------
4972 Stuff to read a file in "delay" mode -- 01 Jan 1997.
4973 -----------------------------------------------------------------------------*/
4974
4975 /**** If possible, throw the data away for later retrieval from disk ****/
4976
mri_purge_delay(MRI_IMAGE * im)4977 void mri_purge_delay( MRI_IMAGE * im )
4978 {
4979 void * ar ;
4980
4981 /** if no delay filename,
4982 or if it is marked as already set for delay input, do nothing **/
4983
4984 if( im->fname == NULL ||
4985 (im->fondisk & INPUT_DELAY) != 0 ) return ;
4986
4987 /** get the data pointer, throw data way, clear the data pointer **/
4988
4989 ar = mri_data_pointer( im ) ;
4990 if( ar != NULL ){ free(ar) ; mri_clear_data_pointer(im) ; }
4991
4992 /** mark as set for delay input **/
4993
4994 im->fondisk |= INPUT_DELAY ;
4995 return ;
4996 }
4997
4998 /**** if possible, read data from delay input file ****/
4999
mri_input_delay(MRI_IMAGE * im)5000 void mri_input_delay( MRI_IMAGE * im )
5001 {
5002 FILE * imfile=NULL ;
5003 void * imar ;
5004
5005 /** if no delay input file,
5006 or is marked as already read in, do nothing **/
5007
5008 if( im->fname == NULL ||
5009 (im->fondisk & INPUT_DELAY) == 0 ) return ;
5010
5011 /** open the delay input file [06 Mar 2001: maybe not] **/
5012
5013 if( strcmp(im->fname,"ALLZERO") != 0 ){
5014 imfile = fopen( im->fname , "r" ) ;
5015 if( imfile == NULL ){
5016 ERROR_message("mri_input_delay: couldn't open file %s" , im->fname ) ; return ;
5017 }
5018 }
5019
5020 /** make space for the array **/
5021
5022 imar = (void *) malloc( im->nvox * im->pixel_size ) ;
5023 if( imar == NULL ){
5024 fprintf( stderr ,
5025 "malloc fails for delayed image from file %s\n" , im->fname ) ;
5026 if( imfile != NULL ) fclose( imfile ) ;
5027 return ;
5028 }
5029 mri_fix_data_pointer( imar , im ) ;
5030
5031 /** read from the file into the array **/
5032
5033 if( imfile != NULL ){
5034 fseek( imfile , im->foffset , SEEK_SET ) ;
5035 fread( imar , im->pixel_size , im->nvox , imfile ) ;
5036 fclose( imfile ) ;
5037 } else {
5038 memset( imar , 0 , im->nvox * im->pixel_size ) ; /* 06 Mar 2001 */
5039 }
5040
5041 /** swap bytes, if so marked **/
5042
5043 if( (im->fondisk & BSWAP_DELAY) ){
5044 mri_swapbytes( im ) ;
5045 im->was_swapped = 1 ; /* 07 Mar 2002 */
5046 }
5047
5048 /** mark as already read from disk **/
5049
5050 im->fondisk ^= INPUT_DELAY ;
5051
5052 #if 0
5053 fprintf(stderr,"delayed input from file %s at offset %d\n",im->fname,im->foffset);
5054 #endif
5055 return ;
5056 }
5057
5058 /**********************************************************************/
5059 /**** like mri_read_file, but returns delayed images for 3D: files ****/
5060 /**** (all others are read in now anyhoo, so there) ****/
5061
mri_read_file_delay(char * fname)5062 MRI_IMARR * mri_read_file_delay( char * fname )
5063 {
5064 MRI_IMARR *newar=NULL ;
5065 MRI_IMAGE *newim ;
5066 char *new_fname ;
5067 int tried_dicom=0 ;
5068
5069 new_fname = imsized_fname( fname ) ;
5070 if( new_fname == NULL ) return NULL ;
5071
5072 if( strlen(new_fname) > 9 && new_fname[0] == '3' && new_fname[1] == 'D' &&
5073 (new_fname[2] == ':' || new_fname[3] == ':') ){
5074 /* check for ':', too 3 Jan 2005 [rickr] */
5075
5076 newar = mri_read_3D_delay( new_fname ) ; /* read from a 3D file, later */
5077
5078 } else if( strlen(new_fname) > 9 &&
5079 new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
5080
5081 newar = mri_read_3A( new_fname ) ;
5082
5083 } else if( check_dicom_magic_num( new_fname ) ) {
5084
5085 newar = mri_read_dicom( new_fname ); tried_dicom=1 ;
5086
5087 } else if( strstr(new_fname,".hdr") != NULL ||
5088 strstr(new_fname,".HDR") != NULL ){ /* 05 Feb 2001 - ANALYZE header */
5089
5090 newar = mri_read_analyze75( new_fname ) ;
5091
5092 } else if( strstr(new_fname,".ima") != NULL ||
5093 strstr(new_fname,".IMA") != NULL ){ /* 12 Mar 2001 - Siemens */
5094
5095 newar = mri_read_siemens( new_fname ) ;
5096
5097 } else if( strstr(new_fname,".mpg" ) != NULL || /* 03 Dec 2003 */
5098 strstr(new_fname,".MPG" ) != NULL || /* read MPEGs */
5099 strstr(new_fname,".mpeg") != NULL ||
5100 strstr(new_fname,".MPEG") != NULL ){
5101
5102 newar = mri_read_mpeg( new_fname ) ; /* cf. mri_read_mpeg.c */
5103 }
5104
5105 /* failed thus far? try DICOM, unless user has requested DICOM last */
5106 /* 05 May 2003 added option to try DICOM last KRH */
5107
5108 if ((newar == NULL) && !AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
5109 if( !tried_dicom ){
5110 newar = mri_read_dicom( new_fname ) ; tried_dicom = 1 ;
5111 }
5112 }
5113
5114 /* failed again? try mri_read() for 1 image */
5115
5116 if( newar == NULL ){
5117 newim = mri_read( new_fname ) ; /* read from a 2D file */
5118 if( newim == NULL ){ free(new_fname) ; return NULL ; }
5119 INIT_IMARR(newar) ;
5120 ADDTO_IMARR(newar,newim) ;
5121 }
5122
5123 if ( (newar == NULL) && AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
5124 if( !tried_dicom ){
5125 newar = mri_read_dicom( new_fname ) ; tried_dicom = 1 ;
5126 }
5127 }
5128
5129 free(new_fname) ;
5130 return newar ;
5131 }
5132
5133 /**** like mri_read_3D, but returns delayed images ****/
5134
mri_read_3D_delay(char * tname)5135 MRI_IMARR * mri_read_3D_delay( char * tname )
5136 {
5137 int hglobal , himage , nx , ny , nz ;
5138 char fname[256] , buf[512] ;
5139 int ngood , kim , datum_type , datum_len , swap ;
5140 MRI_IMARR *newar ;
5141 MRI_IMAGE *newim ;
5142 FILE *imfile ;
5143 long long length , nneed , hglob=0 ; /* 22 Mar 2007 */
5144
5145 /*** get info from 3D tname ***/
5146
5147 if( tname == NULL || strlen(tname) < 10 ) return NULL ;
5148
5149 switch( tname[2] ){ /* allow for 3D: or 3Ds: or 3Db: */
5150
5151 default:
5152 case ':':
5153 ngood = sscanf( tname , "3D:%d:%d:%d:%d:%d:%s" ,
5154 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5155
5156 swap = 0 ;
5157 datum_type = MRI_short ;
5158 datum_len = sizeof(short) ; /* better be 2 */
5159 break ;
5160
5161 case 's':
5162 ngood = sscanf( tname , "3Ds:%d:%d:%d:%d:%d:%s" ,
5163 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5164
5165 swap = 1 ;
5166 datum_type = MRI_short ;
5167 datum_len = sizeof(short) ; /* better be 2 */
5168 break ;
5169
5170 case 'b':
5171 ngood = sscanf( tname , "3Db:%d:%d:%d:%d:%d:%s" ,
5172 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5173
5174 swap = 0 ;
5175 datum_type = MRI_byte ;
5176 datum_len = sizeof(byte) ; /* better be 1 */
5177 break ;
5178
5179 case 'f':
5180 ngood = sscanf( tname , "3Df:%d:%d:%d:%d:%d:%s" ,
5181 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5182
5183 swap = 0 ;
5184 datum_type = MRI_float ;
5185 datum_len = sizeof(float) ; /* better be 4 */
5186 break ;
5187
5188 case 'd': /* 06 Feb 2003 */
5189 ngood = sscanf( tname , "3Dd:%d:%d:%d:%d:%d:%s" ,
5190 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5191
5192 swap = 0 ;
5193 datum_type = MRI_float ;
5194 datum_len = sizeof(double) ; /* better be 8 */
5195 break ;
5196
5197 case 'i':
5198 ngood = sscanf( tname , "3Di:%d:%d:%d:%d:%d:%s" ,
5199 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5200
5201 swap = 0 ;
5202 datum_type = MRI_int ;
5203 datum_len = sizeof(int) ; /* better be 4 */
5204 break ;
5205
5206 case 'c':
5207 ngood = sscanf( tname , "3Dc:%d:%d:%d:%d:%d:%s" ,
5208 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5209
5210 swap = 0 ;
5211 datum_type = MRI_complex ;
5212 datum_len = sizeof(complex) ; /* better be 8 */
5213 break ;
5214
5215 case 'r':
5216 ngood = sscanf( tname , "3Dr:%d:%d:%d:%d:%d:%s" ,
5217 &hglobal , &himage , &nx , &ny , &nz , fname ) ;
5218
5219 swap = 0 ;
5220 datum_type = MRI_rgb ;
5221 datum_len = 3*sizeof(byte) ; /* better be 3 */
5222 break ;
5223 }
5224
5225 if( ngood < 6 || himage < 0 ||
5226 nx <= 0 || ny <= 0 || nz <= 0 ||
5227 strlen(fname) <= 0 ) return NULL ; /* bad info */
5228
5229 /*** open the input file and position it [06 Mar 2001: maybe not] ***/
5230
5231 if( strcmp(fname,"ALLZERO") != 0 ){
5232 imfile = fopen( fname , "r" ) ;
5233 if( imfile == NULL ){
5234 ERROR_message("mri_read_3D_delay: couldn't open file %s" , fname ) ; return NULL ;
5235 }
5236 } else {
5237 imfile = NULL ;
5238 }
5239
5240 if( imfile != NULL ){
5241 length = THD_filesize(fname) ; /* 22 Mar 2007 */
5242
5243 /** 13 Apr 1999: modified to allow actual hglobal < -1
5244 as long as hglobal+himage >= 0 **/
5245
5246 hglob = hglobal ;
5247 if( hglob == -1 || hglob+himage < 0 ){
5248 hglob = length - nz*(datum_len*nx*ny+himage) ;
5249 if( hglob < 0 ) hglob = 0 ;
5250 }
5251
5252 nneed = hglob + (datum_len*nx*ny+himage) * (long long)nz ;
5253 if( length < nneed ){
5254 fprintf( stderr,
5255 "file %s is %lld bytes long but must be at least %lld bytes long\n"
5256 "for hglobal=%lld himage=%d nx=%d ny=%d nz=%d and voxel=%d bytes\n",
5257 fname,length,nneed,hglob,himage,nx,ny,nz,datum_len ) ;
5258 fclose( imfile ) ;
5259 return NULL ;
5260 }
5261 fclose( imfile ) ;
5262 }
5263
5264 /*** put pointers to data in the file into the images ***/
5265
5266 INIT_IMARR(newar) ;
5267
5268 for( kim=0 ; kim < nz ; kim++ ){
5269 newim = mri_new_vol_empty( nx,ny,1 , datum_type ) ; /* empty image */
5270 mri_add_fname_delay( fname , newim ) ; /* put filename in */
5271 newim->fondisk = (swap) ? (INPUT_DELAY | BSWAP_DELAY) /* mark read type */
5272 : (INPUT_DELAY) ;
5273 newim->foffset = hglob + (kim+1)*himage + datum_len*nx*ny*(long long)kim ;
5274
5275 if( nz == 1 ) mri_add_name( fname , newim ) ;
5276 else {
5277 sprintf( buf , "%s#%d" , fname,kim ) ;
5278 mri_add_name( buf , newim ) ;
5279 }
5280
5281 ADDTO_IMARR(newar,newim) ;
5282 }
5283
5284 return newar ;
5285 }
5286
5287 /*----------------------------------------------------------------------------*/
5288 /* Read the first non-trivial line from a file and parse out the
5289 sub-strings (for use as column headers later), allowing for
5290 column indexes [0,3,$] -- 18 May 2010 [RWCox]
5291 *//*--------------------------------------------------------------------------*/
5292
mri_read_1D_headerline(char * fname)5293 THD_string_array * mri_read_1D_headerline( char *fname )
5294 {
5295 THD_string_array *sar ;
5296 char *dname , *cpt , *dpt , lbuf[NLL] ;
5297 int ii , ibot , nlab , *ivlist ;
5298 FILE *fts ;
5299
5300 ENTRY("mri_read_1D_headerline") ;
5301
5302 /* check for bad-ositiness */
5303
5304 if( fname == NULL || *fname == '\0' ) RETURN(NULL) ;
5305 ii = strlen(fname) ;
5306 if( (ii <= 2 && fname[0] == '-') ||
5307 (ii <= 6 && strncmp(fname,"stdin" ,5) == 0) ||
5308 (ii <= 9 && strncmp(fname,"/dev/fd0",8) == 0) ) RETURN(NULL) ;
5309 if( strncmp(fname,"1D:",3) == 0 ) RETURN(NULL) ;
5310
5311 /* copy filename to local variable for editing purposes */
5312
5313 dname = strdup(fname) ; if( dname[ii-1] == '\'' ) dname[ii-1] = '\0' ;
5314
5315 /*-- split filename and subvector list --*/
5316
5317 cpt = strstr(fname,"[") ; /* column list */
5318 dpt = strstr(fname,"{") ; /* we don't use this row list */
5319
5320 if( cpt == fname || dpt == fname ){ /* can't be at start of filename! */
5321 free(dname) ; RETURN(NULL) ;
5322 } else { /* got a subvector list */
5323 if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; }
5324 if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; }
5325 }
5326
5327 /* open input file */
5328
5329 fts = fopen(dname,"r") ;
5330 if( fts == NULL ){
5331 ERROR_message("mri_read_1D_headerline: couldn't open file %s" , dname ) ;
5332 free(dname) ; RETURN(NULL) ;
5333 }
5334 free(dname) ;
5335
5336 /* read lines until we get a useful line */
5337
5338 while(1){
5339 lbuf[0] = '\0' ;
5340 dpt = afni_fgets( lbuf , NLL , fts ) ; /* read a line of data */
5341 if( dpt == NULL ){ fclose(fts); RETURN(NULL); } /* error */
5342 ii = strlen(lbuf) ; if( ii == 0 ) continue ; /* nada => loopback */
5343 if( ii == 1 && isspace(lbuf[0]) ) continue ; /* 1 blank only => loopback */
5344 ibot = (lbuf[0] == '#') ? 1 : 0 ; /* start of scanning */
5345 for( ii=ibot ; isspace(lbuf[ii]) ; ii++ ) continue ; /* skip blanks */
5346 if( lbuf[ii] == '\0' ) continue ; /* all blanks => loopback */
5347 ibot = ii ; break ; /* can process this line */
5348 }
5349 fclose(fts) ; /* don't need file, no more, no how, no way */
5350
5351 /* scan line, break into sub-strings (at whitespace), store sub-strings */
5352
5353 INIT_SARR(sar) ; /* output array of strings */
5354 while(1){
5355 for( ii=ibot ; lbuf[ii] != '\0' && !isspace(lbuf[ii]) ; ii++ ) ; /*nada*/
5356 if( ii == ibot ) break ; /* done */
5357 dpt = (char *)malloc(sizeof(char)*(ii-ibot+2)) ; /* save result */
5358 memcpy( dpt , lbuf+ibot , ii-ibot ) ; dpt[ii-ibot] = '\0' ;
5359 ADDTO_SARR(sar,dpt) ;
5360 if( lbuf[ii] == '\0' ) break ; /* end of input line => done */
5361 for( ; isspace(lbuf[ii]) ; ii++ ) ; /* skip whitespace in input line */
5362 if( lbuf[ii] == '\0' ) break ; /* at end of string => done */
5363 ibot = ii ; continue ; /* otherwise, loopback for next sub-string */
5364 }
5365
5366 /* didn't get nuthin? */
5367
5368 nlab = SARR_NUM(sar) ; if( nlab <= 0 ){ DESTROY_SARR(sar); RETURN(NULL); }
5369
5370 /* if have a [sub-column] chooser, use it now to select only some strings */
5371
5372 ivlist = MCW_get_intlist( nlab , cpt ) ;
5373 if( ivlist != NULL && ivlist[0] > 0 ){
5374 int nts=ivlist[0] , *ivl=ivlist+1 ; THD_string_array *qar ;
5375 INIT_SARR(qar) ;
5376 for( ii=0 ; ii < nts ; ii++ ){
5377 if( ivl[ii] >= 0 && ivl[ii] < nlab ){
5378 dpt = strdup( SARR_STRING(sar,ivl[ii]) ) ;
5379 ADDTO_SARR(qar,dpt) ;
5380 }
5381 }
5382 free(ivlist) ; DESTROY_SARR(sar) ; sar = qar ;
5383 }
5384
5385 RETURN(sar) ; /* et viola */
5386 }
5387