00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*! \file
00008     This file contains all the functions for reading image files.
00009     It is primarily used by to3d.c.  It reads 2D images (into MRI_IMAGE struct)
00010     and arrays of 2D images (into MRI_IMARR struct).
00011 */
00012 
00013 #include <stdio.h>
00014 #include <ctype.h>
00015 #include <string.h>
00016 #include <sys/stat.h>
00017 
00018 /*** for non ANSI compilers ***/
00019 
00020 #ifndef SEEK_END
00021 #define SEEK_END 2
00022 #endif
00023 
00024 #ifndef SEEK_SET
00025 #define SEEK_SET 0
00026 #endif
00027 
00028 #include "mrilib.h"
00029 #include "ge4_header.h"
00030 
00031 /* prototype for local function */
00032 short check_dicom_magic_num( char * );
00033 
00034 /*! Global variable to signal image orientation, if possible. */
00035 
00036 char MRILIB_orients[8] = "\0\0\0\0\0\0\0\0" ;  /* 12 Mar 2001 */
00037 
00038 /*! Global variable to signal image slice offset, if possible. */
00039 
00040 float MRILIB_zoff      = 0.0 ;
00041 
00042 /*! Global variable to signal image TR, if possible. */
00043 
00044 float MRILIB_tr        = 0.0 ;   /* 03 Dec 2001 */
00045 
00046 /*! Global variable to signal image x offset, if possible. */
00047 
00048 float MRILIB_xoff      = 0.0 ;   /* 07 Dec 2001 */
00049 
00050 /*! Global variable to signal image y offset, if possible. */
00051 
00052 float MRILIB_yoff      = 0.0 ;
00053 
00054 /*! Global variable saying whether to use MRILIB_xoff. */
00055 
00056 int use_MRILIB_xoff    = 0 ;
00057 
00058 /*! Global variable saying whether to use MRILIB_yoff. */
00059 
00060 int use_MRILIB_yoff    = 0 ;
00061 
00062 /*! Global variable saying whether to use MRILIB_zoff. */
00063 
00064 int use_MRILIB_zoff    = 0 ;
00065 
00066 /*! Global variable saying whether to use MRILIB_xcos. */
00067 
00068 int use_MRILIB_xcos    = 0 ;
00069 
00070 /*! Global vector pointing in direction of x-axis. */
00071 
00072 float MRILIB_xcos[3]   = { 1.0 , 0.0 , 0.0 } ;
00073 
00074 /*! Global variable saying whether to use MRILIB_ycos. */
00075 
00076 int use_MRILIB_ycos    = 0 ;
00077 
00078 /*! Global vector pointing in direction of y-axis. */
00079 
00080 float MRILIB_ycos[3]   = { 0.0 , 1.0 , 0.0 } ;
00081 
00082 /*! Global variable saying whether to use MRILIB_zcos. */
00083 
00084 int use_MRILIB_zcos    = 0 ;
00085 
00086 /*! Global vector pointing in direction of y-axis. */
00087 
00088 float MRILIB_zcos[3]   = { 0.0 , 0.0 , 1.0 } ;
00089 
00090 /*! Global variable saying whether to use MRILIB_slicespacing. */
00091 
00092 int use_MRILIB_slicespacing = 0 ;
00093 
00094 /*! Global variable giving the spacing between slice centers. */
00095 
00096 float MRILIB_slicespacing = 0.0 ;
00097 
00098 /*** 7D SAFE (but most routines only return 2D images!) ***/
00099 
00100 MRI_IMAGE *mri_try_mri( FILE * , int * ) ;  /* prototypes */
00101 MRI_IMAGE *mri_try_7D ( FILE * , int * ) ;
00102 MRI_IMAGE *mri_try_pgm( FILE * , int * ) ;
00103 
00104 /*-----------------------------------------------------------------*/
00105 
00106 /*! Database of preset file sizes for auto-use of 3D.
00107 
00108   If( head < 0 ) then file length must match size
00109   else file length must == n*size + head for some
00110        integer n, which will be the number of slices
00111        to read from the file.
00112 */
00113 
00114 typedef struct {
00115    int size ;       /*!< file size in bytes if head < 0; image size in bytes if head >= 0 */
00116    int head ;       /*!< global header size if >= 0 */
00117    char * prefix ;  /*!< character string to prefix to filename */
00118 } MCW_imsize ;
00119 
00120 /*! Max number of preset file sizes to allow. */
00121 
00122 #define MAX_MCW_IMSIZE 99
00123 
00124 /*! Array of preset file sizes to use when reading image files. */
00125 
00126 static MCW_imsize imsize[MAX_MCW_IMSIZE] ;
00127 
00128 /*! If this < 0 ==> must initialize array of preset file sizes. */
00129 
00130 static int MCW_imsize_good = -1 ;
00131 
00132 /*---------------------------------------------------------------*/
00133 
00134 #undef swap_4
00135 #undef swap_2
00136 
00137 /*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
00138 
00139 static void swap_4(void *ppp)
00140 {
00141    unsigned char *pntr = (unsigned char *) ppp ;
00142    unsigned char b0, b1, b2, b3;
00143 
00144    b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00145    *pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
00146 }
00147 
00148 /*---------------------------------------------------------------*/
00149 
00150 /*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
00151 
00152 static void swap_8(void *ppp)
00153 {
00154    unsigned char *pntr = (unsigned char *) ppp ;
00155    unsigned char b0, b1, b2, b3;
00156    unsigned char b4, b5, b6, b7;
00157 
00158    b0 = *pntr    ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00159    b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
00160 
00161    *pntr     = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
00162    *(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
00163 }
00164 
00165 /*---------------------------------------------------------------*/
00166 
00167 /*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
00168 
00169 static void swap_2(void *ppp)
00170 {
00171    unsigned char *pntr = (unsigned char *) ppp ;
00172    unsigned char b0, b1;
00173 
00174    b0 = *pntr; b1 = *(pntr+1);
00175    *pntr = b1; *(pntr+1) = b0;
00176 }
00177 
00178 /******************************************************************/
00179 
00180 /*! Earliest image reading function in the AFNI package.
00181     Reads a single 2D image.
00182 
00183     \param  fname is the name of the file to try to read
00184     \return NULL if an image couldn't be read, otherwise
00185             a pointer to an MRI_IMAGE with data, dimensions, etc.
00186 */
00187 
00188 MRI_IMAGE *mri_read( char *fname )
00189 {
00190    FILE      *imfile ;
00191    MRI_IMAGE *im ;
00192    int       length , skip=0 , swap=0 ;
00193    void      *data ;
00194 
00195 ENTRY("mri_read") ;
00196 
00197    if( fname == NULL || *fname == '\0' ) RETURN(NULL) ;  /* bad user */
00198 
00199    /**-- 27 Apr 2005: check here for special filenames --**/
00200 
00201    if( strstr(fname,".jpg" ) != NULL ||  /* various formats  */
00202        strstr(fname,".JPG" ) != NULL ||  /* that we convert  */
00203        strstr(fname,".jpeg") != NULL ||  /* to PPG/PGM using */
00204        strstr(fname,".JPEG") != NULL ||  /* external filters */
00205        strstr(fname,".gif" ) != NULL ||
00206        strstr(fname,".GIF" ) != NULL ||
00207        strstr(fname,".tif" ) != NULL ||
00208        strstr(fname,".TIF" ) != NULL ||
00209        strstr(fname,".tiff") != NULL ||
00210        strstr(fname,".TIFF") != NULL ||
00211        strstr(fname,".bmp" ) != NULL ||
00212        strstr(fname,".BMP" ) != NULL ||
00213        strstr(fname,".pbm" ) != NULL ||
00214        strstr(fname,".PBM" ) != NULL ||
00215        strstr(fname,".pgm" ) != NULL ||
00216        strstr(fname,".PGM" ) != NULL ||
00217        strstr(fname,".ppm" ) != NULL ||
00218        strstr(fname,".PPM" ) != NULL ||
00219        strstr(fname,".png" ) != NULL ||
00220        strstr(fname,".PNG" ) != NULL   ){
00221 
00222      im = mri_read_stuff(fname) ; if( im != NULL ) RETURN(im) ;
00223    }
00224 
00225    /*-- check if file exists and is readable --*/
00226 
00227    imfile = fopen( fname , "r" ) ;
00228    if( imfile == NULL ){
00229      fprintf( stderr , "couldn't open image file %s\n" , fname ) ;
00230      RETURN( NULL );
00231    }
00232 
00233    fseek( imfile , 0L , SEEK_END ) ;  /* get the length of the file */
00234    length = ftell( imfile ) ;         /* (the AJ way) */
00235 
00236    /*--- 03 Dec 2001: check for GEMS format file "IMGF"   ---*/
00237    /*[[[ Information herein from Medical Image Format FAQ ]]]*/
00238 
00239    { char str[5]="AFNI" ;
00240      int nx , ny , bpp , cflag , hdroff , extraskip=0 ;
00241      rewind(imfile) ; fread(str,1,4,imfile) ;   /* check for "IMGF" or "GEMS" */
00242 
00243      if( str[0]=='G' && str[1]=='E' && str[2]=='M' && str[3]=='S' ){ /* 12 Feb 2004 */
00244        char buf[4096]; int bb,cc;                /* search for IMGF in 1st 4K */
00245        rewind(imfile); cc=fread(buf,1,4096,imfile); cc-=4 ;
00246        for( bb=4; bb < cc ; bb++ )
00247         if( buf[bb]=='I' && buf[bb+1]=='M' && buf[bb+2]=='G' && buf[bb+3]=='F' ) break ;
00248        if( bb < cc ){
00249          fseek( imfile , (long)bb , SEEK_SET ) ; extraskip = bb ;
00250          fread(str,1,4,imfile) ;
00251        }
00252      }
00253 
00254      /* 12 Feb 2004: modified to allow for starting at extraskip */
00255 
00256      if( str[0]=='I' && str[1]=='M' && str[2]=='G' && str[3]=='F' ){
00257 
00258        fread( &skip , 4,1, imfile ) ;  /* read next 5 ints */
00259        fread( &nx   , 4,1, imfile ) ;
00260        fread( &ny   , 4,1, imfile ) ;
00261        fread( &bpp  , 4,1, imfile ) ;
00262        fread( &cflag, 4,1, imfile ) ;
00263 
00264        if( nx < 0 || nx > 8192 ){      /* maybe have to byte swap 5 ints */
00265          swap = 1 ;                    /* flag that we are swapping data */
00266          swap_4(&skip); swap_4(&nx) ;
00267          swap_4(&ny)  ; swap_4(&bpp); swap_4(&cflag);
00268        }
00269        skip += extraskip ;             /* location of image data in file */
00270        if( nx < 0 || nx > 8192 || ny < 0 || ny > 8192 ) goto The_Old_Way ;
00271        if( skip < 0  || skip  >= length )               goto The_Old_Way ;
00272        if( bpp != 16 || cflag != 1      )               goto The_Old_Way ;
00273 
00274        /* make image space */
00275 
00276        im = mri_new( nx , ny , MRI_short ) ;
00277 
00278        /* try to read image auxiliary data as well (not mandatory) */
00279 
00280        length = fseek( imfile , 148L+extraskip , SEEK_SET ) ; /* magic GEMS offset */
00281        if( length == 0 ){
00282           fread( &hdroff , 4,1 , imfile ) ;  /* location of image header */
00283           if( swap ) swap_4(&hdroff) ;
00284           if( hdroff > 0 ){                  /* read from image header */
00285              float dx,dy,dz, dxx,dyy,dzz, xyz[9], zz ; int itr, ii,jj,kk, qq ;
00286              static int nzoff=0 ;
00287              static float zoff ;
00288 
00289              /* get voxel grid sizes */
00290 
00291              fseek( imfile , hdroff+26+extraskip , SEEK_SET ) ;
00292              fread( &dzz , 4,1 , imfile ) ;
00293 
00294              fseek( imfile , hdroff+50+extraskip , SEEK_SET ) ;
00295              fread( &dxx , 4,1 , imfile ) ;
00296              fread( &dyy , 4,1 , imfile ) ;
00297 
00298              if( swap ){ swap_4(&dxx); swap_4(&dyy); swap_4(&dzz); }
00299 
00300              /* save into image header [dw > 0 is signal that dx,dy,dz are OK] */
00301 
00302              if( dxx > 0.01 && dyy > 0.01 && dzz > 0.01 ){
00303                im->dx = dxx; im->dy = dyy; im->dz = dzz; im->dw = 1.0;
00304              }
00305 
00306              /* grid orientation: from 3 sets of LPI corner coordinates: */
00307              /*   xyz[0..2] = top left hand corner of image     (TLHC)   */
00308              /*   xyz[3..5] = top right hand corner of image    (TRHC)   */
00309              /*   xyz[6..8] = bottom right hand corner of image (BRHC)   */
00310              /* GEMS coordinate orientation here is LPI.                 */
00311              /* Orientation is saved into global string MRILIB_orients.  */
00312              /* N.B.: AFNI coordinates are RAI orientation.              */
00313 
00314              fseek( imfile , hdroff+154+extraskip , SEEK_SET ) ;
00315              fread( xyz , 4,9 , imfile ) ;
00316              if( swap ) swap_fourbytes(9,xyz) ;
00317 
00318              /* x-axis orientation */
00319              /* ii determines which spatial direction is x-axis  */
00320              /* and is the direction that has the biggest change */
00321              /* between the TLHC and TRHC                        */
00322 
00323              dx = fabs(xyz[3]-xyz[0]) ; ii = 1 ;
00324              dy = fabs(xyz[4]-xyz[1]) ; if( dy > dx ){ ii=2; dx=dy; }
00325              dz = fabs(xyz[5]-xyz[2]) ; if( dz > dx ){ ii=3;        }
00326              dx = xyz[ii+2]-xyz[ii-1] ; if( dx < 0. ){ ii = -ii;    }
00327              switch( ii ){
00328                case  1: MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; break;
00329                case -1: MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; break;
00330                case  2: MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; break;
00331                case -2: MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; break;
00332                case  3: MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; break;
00333                case -3: MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; break;
00334                default: MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; break;
00335              }
00336 
00337              /* y-axis orientation */
00338              /* jj determines which spatial direction is y-axis  */
00339              /* and is the direction that has the biggest change */
00340              /* between the BRHC and TRHC                        */
00341 
00342              dx = fabs(xyz[6]-xyz[3]) ; jj = 1 ;
00343              dy = fabs(xyz[7]-xyz[4]) ; if( dy > dx ){ jj=2; dx=dy; }
00344              dz = fabs(xyz[8]-xyz[5]) ; if( dz > dx ){ jj=3;        }
00345              dx = xyz[jj+5]-xyz[jj+2] ; if( dx < 0. ){ jj = -jj;    }
00346              switch( jj ){
00347                case  1: MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; break;
00348                case -1: MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; break;
00349                case  2: MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; break;
00350                case -2: MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; break;
00351                case  3: MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; break;
00352                case -3: MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; break;
00353                default: MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; break;
00354              }
00355 
00356              MRILIB_orients[6] = '\0' ;   /* terminate orientation string */
00357 
00358              kk = 6 - abs(ii)-abs(jj) ;   /* which spatial direction is z-axis */
00359                                           /* where 1=L-R, 2=P-A, 3=I-S */
00360              zz = xyz[kk-1] ;             /* z-coordinate of this slice (from TLHC) */
00361 
00362              im->zo = zz ;                /* 07 Aug 2002: save slice offset */
00363 
00364              /* getting orientation of z-axis requires 2 images in a row -*/
00365 
00366              if( nzoff == 0 ){  /* from 1st GEMS image */
00367 
00368                zoff = zz ;      /* save this for 2nd image calculation */
00369                switch( kk ){    /* may be changed on second image */
00370                 case  1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
00371                 case  2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
00372                 case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
00373                 default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
00374                }
00375 
00376              } else if( nzoff == 1 ){   /* from 2nd GEMS image */
00377 
00378                float qoff = zz - zoff ;  /* vive la difference */
00379                if( qoff < 0 ) kk = -kk ; /* kk determines z-axis orientation */
00380 
00381                if( !use_MRILIB_slicespacing && qoff != 0.0 ){ /* 10 Jan 2004 */
00382                  use_MRILIB_slicespacing = 1 ;
00383                      MRILIB_slicespacing = fabs(qoff) ;
00384                }
00385 
00386                switch( kk ){
00387                 case  1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
00388                 case -1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
00389                 case  2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
00390                 case -2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
00391                 case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
00392                 case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
00393                 default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
00394                }
00395 
00396                /* save spatial offset of first slice              */
00397                /* [this needs to be positive in the direction of] */
00398                /* [the -z axis, so may need to change its sign  ] */
00399 
00400                MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
00401                if( kk == 1 || kk == 2 || kk == 3 ) MRILIB_zoff = -MRILIB_zoff ;
00402 
00403                /* Same for x offset; [20 Dec 2001]
00404                   This must be at the middle of the TLHC voxel,
00405                     so we must move a little bit towards the TRHC edge;
00406                   We only use the result if the x-coordinate doesn't
00407                     change significantly between the TRHC and BRHC,
00408                     to avoid problems with oblique slices.         */
00409 
00410                qq = abs(ii) ;
00411                MRILIB_xoff = ( xyz[qq-1]*(nx-0.5) + xyz[qq+2]*0.5 ) / nx ;
00412                if( ii == 1 || ii == 2 || ii == 3 ) MRILIB_xoff = -MRILIB_xoff ;
00413                use_MRILIB_xoff = ( fabs(xyz[qq+2]-xyz[qq+5]) < 0.01*dxx ) ;
00414 
00415                /* Same for y offset;
00416                   This must be at the middle of the TLHC voxel,
00417                     so we must move a little bit towards the BRHC edge;
00418                   We only use the result if the y-coordinate doesn't
00419                     change significantly between the TLHC and TRHC. */
00420 
00421                qq = abs(jj) ;
00422                MRILIB_yoff = ( xyz[qq-1]*(ny-0.5) + xyz[qq+5]*0.5 ) / ny ;
00423                if( jj == 1 || jj == 2 || jj == 3 ) MRILIB_yoff = -MRILIB_yoff ;
00424                use_MRILIB_yoff = ( fabs(xyz[qq-1]-xyz[qq+2]) < 0.01*dyy ) ;
00425              }
00426              nzoff++ ;  /* 3rd and later images don't count for z-orientation */
00427 
00428              /* get TR, save into global variable */
00429 
00430              if( MRILIB_tr <= 0.0 ){
00431                fseek( imfile , hdroff+194+extraskip , SEEK_SET ) ;
00432                fread( &itr , 4,1 , imfile ) ;
00433                if( swap ) swap_4(&itr) ;
00434                MRILIB_tr = im->dt = 1.0e-6 * itr ;
00435              }
00436           }
00437        } /* end of trying to read image header */
00438 
00439        goto Ready_To_Roll ;  /* skip to the reading place */
00440      }
00441    } /* end of GEMS */
00442 
00443    /*--- OK, do it the old way ---*/
00444 
00445 The_Old_Way:
00446 
00447 #if 0
00448    MRILIB_orients[0] = '\0' ; MRILIB_zoff = MRILIB_tr = 0.0 ;  /* 03 Dec 2001 */
00449 #endif
00450 
00451    switch( length ){
00452 
00453       case 512:    /* raw 16x16 short -- RWCox: 06 Dec 2001 */
00454          im = mri_new( 16 , 16 , MRI_short ) ;
00455          break ;
00456 
00457       case 2048:   /* raw 32x32 short -- RWCox: 19 Sep 2000 */
00458          im = mri_new( 32 , 32 , MRI_short ) ;
00459          break ;
00460 
00461       case 4096:   /* raw 64x64 byte -- RWC 3/21/95 */
00462          im = mri_new( 64 , 64 , MRI_byte ) ;
00463          break ;
00464 
00465       case 8192:   /* raw 64x64 short */
00466       case 16096:  /* with Signa 5.x header */
00467          im = mri_new( 64 , 64 , MRI_short ) ;
00468          skip = length - 8192 ;
00469          break ;
00470 
00471 #if 0
00472       case 18432:  /* raw 96x96 short */
00473          im = mri_new( 96 , 96 , MRI_short ) ;
00474          break ;
00475 #endif
00476 
00477       case 16384:  /* raw 128x128 byte -- RWC 3/21/95 */
00478          im = mri_new( 128 , 128 , MRI_byte ) ;
00479          break ;
00480 
00481       case 32768:  /* raw 128x128 short */
00482       case 40672:  /* with Signa 5.x header */
00483          im = mri_new( 128 , 128 , MRI_short ) ;
00484          skip = length - 32768 ;
00485          break ;
00486 
00487       case 65536:  /* raw 256x256 8-bit -- Matthew Belmonte March 1995 */
00488          im = mri_new( 256 , 256 , MRI_byte ) ;
00489          break ;
00490 
00491       case 131072:  /* raw 256x256 short */
00492       case 138976:  /* Signa 5.x */
00493       case 145408:  /* Signa 4.x */
00494 
00495          im   = mri_new( 256 , 256 , MRI_short ) ;
00496          skip = length - 131072 ;
00497          break ;
00498 
00499 #if 0
00500       case 262144:  /* raw 256x256 float */
00501          im = mri_new( 256 , 256 , MRI_float ) ;
00502          break ;
00503 #else
00504       case 262144:  /* raw 512x512 byte -- RWC 3/21/95 */
00505          im = mri_new( 512 , 512 , MRI_byte ) ;
00506          break ;
00507 
00508       case 524288:  /* raw 512x512 short -- RWC 3/21/95 */
00509          im = mri_new( 512 , 512 , MRI_short ) ;
00510          break ;
00511 
00512       case 1048576: /* raw 1024x1024 byte -- RWC 3/21/95 */
00513          im = mri_new( 1024 , 1024 , MRI_byte ) ;
00514          break ;
00515 
00516       case 2097152: /* raw 1024x1024 short -- RWC 3/21/95 */
00517          im = mri_new( 1024 , 1024 , MRI_short ) ;
00518          break ;
00519 #endif
00520 
00521       /** not a canonical length: try something else **/
00522 
00523       default:
00524                           im = mri_try_mri( imfile , &skip ) ;  /* Cox format */
00525          if( im == NULL ) im = mri_try_7D ( imfile , &skip ) ;  /* 7D format  */
00526          if( im == NULL ) im = mri_try_pgm( imfile , &skip ) ;  /* PGM format */
00527          if( im != NULL ) break ;
00528 
00529          fclose( imfile ) ; /* close it, since we failed (so far) */
00530 
00531          im = mri_read_ascii( fname ) ;    /* list of ASCII numbers */
00532          if( im != NULL ) RETURN( im );
00533 
00534          im = mri_read_ppm( fname ) ;      /* 15 Apr 1999 */
00535          if( im != NULL ) RETURN( im );
00536 
00537          im = mri_read_stuff( fname ) ;    /* 22 Nov 2002 */
00538          if( im != NULL ) RETURN( im );
00539 
00540          fprintf( stderr , "do not recognize image file %s\n" , fname );
00541          fprintf( stderr , "length seen as %d\n" , length ) ;
00542          RETURN( NULL );
00543    }
00544 
00545    /*-- Actually read the data from disk --*/
00546 
00547 Ready_To_Roll:
00548 
00549    data = mri_data_pointer( im ) ;
00550 
00551    length = fseek( imfile , skip , SEEK_SET ) ;
00552    if( length != 0 ){
00553       fprintf( stderr , "mri_read error in skipping in file %s\n" , fname ) ;
00554       mri_free( im ) ;
00555       RETURN( NULL );
00556    }
00557 
00558    length = fread( data , im->pixel_size , im->nvox , imfile ) ;
00559    fclose( imfile ) ;
00560 
00561    if( length != im->nvox ){
00562       mri_free( im ) ;
00563       fprintf( stderr , "couldn't read image data from file %s\n" , fname ) ;
00564       RETURN( NULL );
00565    }
00566 
00567    mri_add_name( fname , im ) ;
00568 
00569    /*-- 03 Dec 2001: maybe need to swap bytes --*/
00570 
00571    if( swap ){
00572      switch( im->pixel_size ){
00573        default: break ;
00574        case 2:  swap_twobytes (   im->nvox, data ) ; break ;  /* short */
00575        case 4:  swap_fourbytes(   im->nvox, data ) ; break ;  /* int, float */
00576        case 8:  swap_fourbytes( 2*im->nvox, data ) ; break ;  /* complex */
00577      }
00578 
00579      im->was_swapped = 1 ;  /* 07 Mar 2002 */
00580    }
00581 
00582    RETURN( im );
00583 }
00584 
00585 
00586 /******************************************************************/
00587 
00588 /* GEMS 4.x image reading   - rickr  03 Jun 2003 */
00589 
00590 /*! Read a single 2D GEMS 4.x image.
00591 
00592     \param  filename is the name of the file to try to read
00593     \return NULL if an image could not be read, otherwise
00594             the address of a new MRI_IMAGE structure
00595 */
00596 
00597 MRI_IMAGE * mri_read_ge4( char * filename )
00598 {
00599     MRI_IMAGE * im;
00600     ge4_header  H;
00601 
00602 ENTRY( "mri_read_ge4" );
00603 
00604     if ( filename == NULL )
00605     {
00606         fprintf( stderr, "** mri_read_ge4 - missing filename\n" );
00607         RETURN( NULL );
00608     }
00609 
00610     /* try to read image file - return with image */
00611     if ( ge4_read_header( &H, filename, True ) != 0 )
00612         RETURN( NULL );
00613 
00614     /* these dimensions are fixed */
00615     if ( (im = mri_new(256, 256, MRI_short)) == NULL )
00616     {
00617         free(H.image);
00618         RETURN( NULL );
00619     }
00620 
00621     /* fill im struct with data from H */
00622     im->zo = H.im_h.im_loc;        /* this may well be incorrect */
00623     im->dt = H.im_h.tr;
00624     im->was_swapped = H.swap;
00625 
00626     if ( ( H.ser_h.fov >    1.0      ) &&
00627          ( H.ser_h.fov < 1000.0      ) &&
00628          ( H.ser_h.scan_mat_x >    0 ) &&
00629          ( H.ser_h.scan_mat_x < 1000 ) &&
00630          ( H.ser_h.scan_mat_y >    0 ) &&
00631          ( H.ser_h.scan_mat_y < 1000 ) )
00632     {
00633         /* attempt to set dx, dy and dz from these */
00634 
00635         im->dx = 2 * H.ser_h.fov / H.ser_h.scan_mat_x;
00636         im->dy = im->dx;
00637         im->dz = 2 * H.ser_h.fov / H.ser_h.scan_mat_y;
00638         im->dw = 1;
00639     }
00640 
00641     memcpy( mri_data_pointer(im), H.image, H.im_bytes );
00642 
00643     mri_add_name( filename, im );
00644 
00645     free(H.image);        /* your services are no longer required */
00646 
00647     RETURN( im );
00648 }
00649 
00650 
00651 /*********************************************************************/
00652 
00653 #if 0
00654 #define NUMSCAN(var)                                                       \
00655    { while( isspace(ch) ) {ch = getc(imfile) ;}                            \
00656      if(ch == '#') do{ch = getc(imfile) ;}while(ch != '\n' && ch != EOF) ; \
00657      if(ch =='\n') ch = getc(imfile) ;                                     \
00658      for( nch=0 ; isdigit(ch) ; nch++,ch=getc(imfile) ) {buf[nch] = ch ;}  \
00659      buf[nch]='\0';                                                        \
00660      var = strtol( buf , NULL , 10 ) ; }
00661 #else
00662 
00663 /*! Skip comments in a PPM file. */
00664 
00665 #define SKIPCOM                                                            \
00666     {if(ch == '#') do{ch = getc(imfile) ;}while(ch != '\n' && ch != EOF);}
00667 
00668 /*! Scan for a number in the imfile stream. */
00669 
00670 #define NUMSCAN(var)                                                       \
00671    { SKIPCOM ;                                                             \
00672      while( ch!=EOF && !isdigit(ch) ){ch = getc(imfile); SKIPCOM; }        \
00673      for( nch=0 ; isdigit(ch) ; nch++,ch=getc(imfile) ) {buf[nch] = ch ;}  \
00674      buf[nch]='\0';                                                        \
00675      var = strtol( buf , NULL , 10 ) ; }
00676 #endif
00677 
00678 /*! Try to read an file in the "Cox MRI" format.
00679 
00680     \param imfile is a pointer to an open FILE.
00681     \param skip is a pointer to an int that will be set to the number
00682            of bytes to skip from the file start to find the image data
00683     \return NULL if the file doesn't work for "Cox MRI" format;
00684             otherwise, the return is a pointer to an MRI_IMAGE ready
00685             to have its data read from imfile.
00686 */
00687 
00688 MRI_IMAGE *mri_try_mri( FILE *imfile , int *skip )
00689 {
00690    int ch , nch , nx,ny,imcode ;
00691    char buf[64] ;
00692    MRI_IMAGE *im ;
00693 
00694 ENTRY("mri_try_mri") ;
00695 
00696    fseek( imfile , 0 , SEEK_SET ) ;  /* rewind file */
00697 
00698    ch = getc( imfile ) ; if( ch != 'M' ) RETURN( NULL );  /* check for MRI */
00699    ch = getc( imfile ) ; if( ch != 'R' ) RETURN( NULL );
00700    ch = getc( imfile ) ; if( ch != 'I' ) RETURN( NULL );
00701 
00702    /* magic MRI found, so read numbers */
00703 
00704    ch = getc(imfile) ;
00705 
00706    NUMSCAN(imcode) ;
00707    NUMSCAN(nx) ;  if( nx <= 0 ) RETURN( NULL );
00708    NUMSCAN(ny) ;  if( ny <= 0 ) RETURN( NULL );
00709 
00710    *skip = ftell(imfile) ;
00711    im    = mri_new( nx , ny , imcode ) ;
00712    RETURN( im );
00713 }
00714 
00715 /**************************************************************************
00716    7D format: MRn kind n-dimensions data, where 'n' = 1-7.
00717 ***************************************************************************/
00718 
00719 /*! Try to read a "Cox nD MRI" image file (fat chance). */
00720 
00721 MRI_IMAGE *mri_try_7D( FILE *imfile , int *skip )
00722 {
00723    int ch , nch , nx,ny,nz,nt,nu,nv,nw , imcode , ndim ;
00724    char buf[64] ;
00725    MRI_IMAGE *im ;
00726 
00727 ENTRY("mri_try_7D") ;
00728 
00729    fseek( imfile , 0 , SEEK_SET ) ;  /* rewind file */
00730 
00731    ch = getc( imfile ) ; if( ch != 'M' ) RETURN( NULL );  /* check for MR[1-7] */
00732    ch = getc( imfile ) ; if( ch != 'R' ) RETURN( NULL );
00733    ch = getc( imfile ) ;
00734    switch( ch ){
00735       default:  RETURN( NULL );   /* not what I expected */
00736 
00737       case '1': ndim = 1 ; break ;
00738       case '2': ndim = 2 ; break ;
00739       case '3': ndim = 3 ; break ;
00740       case '4': ndim = 4 ; break ;
00741       case '5': ndim = 5 ; break ;
00742       case '6': ndim = 6 ; break ;
00743       case '7': ndim = 7 ; break ;
00744    }
00745    /* magic MR? found, so read numbers */
00746 
00747    ch = getc(imfile) ;
00748    NUMSCAN(imcode) ;
00749 
00750    nx = ny = nz = nt = nu = nv = nw = 1 ;
00751 
00752                    NUMSCAN(nx) ;  if( nx <= 0 ) RETURN( NULL );
00753    if( ndim > 1 ){ NUMSCAN(ny) ;  if( ny <= 0 ) RETURN( NULL ); }
00754    if( ndim > 2 ){ NUMSCAN(nz) ;  if( nz <= 0 ) RETURN( NULL ); }
00755    if( ndim > 3 ){ NUMSCAN(nt) ;  if( nt <= 0 ) RETURN( NULL ); }
00756    if( ndim > 4 ){ NUMSCAN(nu) ;  if( nu <= 0 ) RETURN( NULL ); }
00757    if( ndim > 5 ){ NUMSCAN(nv) ;  if( nv <= 0 ) RETURN( NULL ); }
00758    if( ndim > 6 ){ NUMSCAN(nw) ;  if( nw <= 0 ) RETURN( NULL ); }
00759 
00760    *skip = ftell(imfile) ;
00761    im    = mri_new_7D_generic( nx,ny,nz,nt,nu,nv,nw , imcode , TRUE ) ;
00762    RETURN( im );
00763 }
00764 
00765 
00766 /*********************************************************************/
00767 
00768 /*! Try to read a raw PGM format image file.
00769 
00770     \param imfile is a pointer to an open FILE
00771     \param skip is a pointer to an int; *skip will be set to the
00772            byte offset at which to start reading data
00773     \return A pointer to an MRI_IMAGE ready to have its data read in
00774             (if the file is a PGM file), or NULL.
00775 */
00776 
00777 MRI_IMAGE *mri_try_pgm( FILE *imfile , int *skip )
00778 {
00779    int ch , nch , nx,ny,maxval ;
00780    char buf[64] ;
00781    MRI_IMAGE *im ;
00782 
00783 ENTRY("mri_try_pgm") ;
00784 
00785    fseek( imfile , 0 , SEEK_SET ) ;  /* rewind file */
00786 
00787    ch = getc( imfile ) ; if( ch != 'P' ) RETURN(NULL);  /* check for magic */
00788    ch = getc( imfile ) ; if( ch != '5' ) RETURN(NULL);
00789 
00790    /* magic P5 found, so read numbers */
00791 
00792    ch = getc(imfile) ;
00793 
00794    NUMSCAN(nx)     ; if( nx     <= 0 ) RETURN(NULL);
00795    NUMSCAN(ny)     ; if( ny     <= 0 ) RETURN(NULL);
00796    NUMSCAN(maxval) ; if( maxval <= 0 || maxval >  255 ) RETURN(NULL);
00797 
00798    *skip = ftell(imfile) ;
00799    im    = mri_new( nx , ny , MRI_byte ) ;
00800    RETURN(im);
00801 }
00802 
00803 /*--------------------------------------------------------------
00804    Read a pile of images from one file.
00805    Modified 4/4/95 to read short or byte data.
00806    Modified 10/02/95 to allow byte swapping with 3Ds:
00807    Modified 11/06/95 to allow float images with 3Df:
00808                  and to allow int images with 3Di:
00809                  and to allow complex images with 3Dc:
00810    Modified 16 Apr 2002 to allow RGB input with 3Dr:
00811 
00812    [N.B.: if this routine is altered, don't forget mri_imcount!]
00813 ----------------------------------------------------------------*/
00814 
00815 /*! Read one or more 2D slices from a "3D:" formatted image file. */
00816 
00817 MRI_IMARR * mri_read_3D( char * tname )
00818 {
00819    int hglobal , himage , nx , ny , nz ;
00820    char fname[256] , buf[512] ;
00821    int ngood , length , kim , koff , datum_type , datum_len , swap ;
00822    MRI_IMARR * newar ;
00823    MRI_IMAGE * newim ;
00824    void      * imar ;
00825    FILE      * imfile ;
00826 
00827 ENTRY("mri_read_3D") ;
00828 
00829    /*** get info from 3D tname ***/
00830 
00831    if( tname == NULL || strlen(tname) < 10 ) RETURN(NULL) ;
00832 
00833    switch( tname[2] ){  /* allow for 3D: or 3Ds: or 3Db:, etc */
00834 
00835       default:
00836       case ':':
00837          ngood = sscanf( tname , "3D:%d:%d:%d:%d:%d:%s" ,
00838                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00839 
00840          swap       = 0 ;
00841          datum_type = MRI_short ;
00842          datum_len  = sizeof(short) ;  /* better be 2 */
00843          break ;
00844 
00845       case 's':
00846          ngood = sscanf( tname , "3Ds:%d:%d:%d:%d:%d:%s" ,
00847                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00848 
00849          swap       = 1 ;
00850          datum_type = MRI_short ;
00851          datum_len  = sizeof(short) ;  /* better be 2 */
00852          break ;
00853 
00854       case 'b':
00855          ngood = sscanf( tname , "3Db:%d:%d:%d:%d:%d:%s" ,
00856                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00857 
00858          swap       = 0 ;
00859          datum_type = MRI_byte ;
00860          datum_len  = sizeof(byte) ;  /* better be 1 */
00861          break ;
00862 
00863       case 'f':
00864          ngood = sscanf( tname , "3Df:%d:%d:%d:%d:%d:%s" ,
00865                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00866 
00867          swap       = 0 ;
00868          datum_type = MRI_float ;
00869          datum_len  = sizeof(float) ;  /* better be 4 */
00870          break ;
00871 
00872       case 'd':                                            /* 06 Feb 2003 */
00873          ngood = sscanf( tname , "3Dd:%d:%d:%d:%d:%d:%s" ,
00874                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00875 
00876          swap       = 0 ;
00877          datum_type = MRI_double ;
00878          datum_len  = sizeof(double) ;  /* better be 8 */
00879          break ;
00880 
00881       case 'i':
00882          ngood = sscanf( tname , "3Di:%d:%d:%d:%d:%d:%s" ,
00883                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00884 
00885          swap       = 0 ;
00886          datum_type = MRI_int ;
00887          datum_len  = sizeof(int) ;  /* better be 4 */
00888          break ;
00889 
00890       case 'c':
00891          ngood = sscanf( tname , "3Dc:%d:%d:%d:%d:%d:%s" ,
00892                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00893 
00894          swap       = 0 ;
00895          datum_type = MRI_complex ;
00896          datum_len  = sizeof(complex) ;  /* better be 8 */
00897          break ;
00898 
00899       case 'r':
00900          ngood = sscanf( tname , "3Dr:%d:%d:%d:%d:%d:%s" ,
00901                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
00902 
00903          swap       = 0 ;
00904          datum_type = MRI_rgb ;
00905          datum_len  = 3*sizeof(byte) ;  /* better be 3 */
00906          break ;
00907    }
00908 
00909    if( ngood < 6 || himage < 0 ||
00910        nx <= 0   || ny <= 0    || nz <= 0 ||
00911        strlen(fname) <= 0                   ) RETURN(NULL);   /* bad info */
00912 
00913    /*** 06 Mar 2001: special case of fname ***/
00914 
00915    if( strcmp(fname,"ALLZERO") == 0 ){
00916       INIT_IMARR(newar) ;
00917       for( kim=0 ; kim < nz ; kim++ ){
00918          newim = mri_new( nx , ny , datum_type ) ;
00919          imar  = mri_data_pointer( newim ) ;
00920          memset( imar , 0 , newim->nvox * newim->pixel_size ) ;
00921          sprintf( buf , "%s#%d" , fname,kim ) ;
00922          mri_add_name( buf , newim ) ;
00923          ADDTO_IMARR(newar,newim) ;
00924       }
00925       RETURN(newar);
00926    }
00927 
00928    /*** open the input file and position it ***/
00929 
00930    imfile = fopen( fname , "r" ) ;
00931    if( imfile == NULL ){
00932       fprintf( stderr , "couldn't open image file %s\n" , fname ) ;
00933       RETURN(NULL);
00934    }
00935 
00936    fseek( imfile , 0L , SEEK_END ) ;  /* get the length of the file */
00937    length = ftell( imfile ) ;
00938 
00939    /** 13 Apr 1999: modified to allow actual hglobal < -1
00940                     as long as hglobal+himage >= 0       **/
00941 
00942 #if 0                 /* the old code */
00943    if( hglobal < 0 ){
00944       hglobal = length - nz*(datum_len*nx*ny+himage) ;
00945       if( hglobal < 0 ) hglobal = 0 ;
00946    }
00947 #else                 /* 13 Apr 1999 */
00948    if( hglobal == -1 || hglobal+himage < 0 ){
00949       hglobal = length - nz*(datum_len*nx*ny+himage) ;
00950       if( hglobal < 0 ) hglobal = 0 ;
00951    }
00952 #endif
00953 
00954    ngood = hglobal + nz*(datum_len*nx*ny+himage) ;
00955    if( length < ngood ){
00956       fprintf( stderr,
00957         "image file %s is %d bytes long but must be at least %d bytes long\n"
00958         "for hglobal=%d himage=%d nx=%d ny=%d nz=%d and voxel=%d bytes\n",
00959         fname,length,ngood,hglobal,himage,nx,ny,nz,datum_len ) ;
00960       fclose( imfile ) ;
00961       RETURN(NULL);
00962    }
00963 
00964    /*** read images from the file ***/
00965 
00966    INIT_IMARR(newar) ;
00967 
00968    for( kim=0 ; kim < nz ; kim++ ){
00969       koff = hglobal + (kim+1)*himage + datum_len*nx*ny*kim ;
00970       fseek( imfile , koff , SEEK_SET ) ;
00971 
00972       newim  = mri_new( nx , ny , datum_type ) ;
00973       imar   = mri_data_pointer( newim ) ;
00974       length = fread( imar , datum_len , nx * ny , imfile ) ;
00975       if( swap ){
00976          mri_swapbytes( newim ) ;
00977          newim->was_swapped = 1 ;  /* 07 Mar 2002 */
00978       }
00979 
00980       if( nz == 1 ) mri_add_name( fname , newim ) ;
00981       else {
00982          sprintf( buf , "%s#%d" , fname,kim ) ;
00983          mri_add_name( buf , newim ) ;
00984       }
00985 
00986       ADDTO_IMARR(newar,newim) ;
00987    }
00988 
00989    fclose(imfile) ;
00990    RETURN(newar);
00991 }
00992 
00993 /*------------------------------------------------------------------------------*/
00994 
00995 /*! Read one or more 2D images from a file.
00996 
00997    This function is the main point of input for to3d.c.
00998    \param fname is the name of the file to read.  This file
00999           might be in one of these formats:
01000            - "3D:" format (implicitly or explicitly)
01001            - "3A:" format
01002            - *.hdr (ANALYZE 2D-4D) format
01003            - *.ima (Siemens 2D array) format
01004            - I.*   (GEMS) format
01005            - PGM format
01006            - PPM format
01007            - GIF, TIFF, JPEG, BMP, PNG formats (thru filters)
01008            - List of ASCII numbers
01009            - pre-defined 2D file size in mri_read()
01010            - "Cox MRI" (if this is what you have, God help you, no one else can)
01011 
01012    \return A pointer to an array of 2D images.  If nothing
01013            could be read, NULL is returned.
01014 */
01015 
01016 MRI_IMARR * mri_read_file( char * fname )
01017 {
01018    MRI_IMARR * newar=NULL ;
01019    MRI_IMAGE * newim ;
01020    char * new_fname ;
01021 
01022 ENTRY("mri_read_file") ;
01023 
01024    /* convert fname to new_fname, based on environment */
01025 
01026    new_fname = imsized_fname( fname ) ;
01027    if( new_fname == NULL ) RETURN( NULL );
01028 
01029    /* input method is based on filename */
01030 
01031    if( strlen(new_fname) > 9 &&
01032        new_fname[0] == '3'   &&
01033        new_fname[1] == 'D'   &&
01034       (new_fname[2] == ':' || new_fname[3] == ':') ){
01035 
01036       newar = mri_read_3D( new_fname ) ;   /* read from a 3D: file */
01037 
01038    } else if( strlen(new_fname) > 9 &&
01039               new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
01040 
01041       newar = mri_read_3A( new_fname ) ;   /* from a 3A: file */
01042 
01043    } else if( check_dicom_magic_num( new_fname ) ) { /* 10 Aug 2004 */
01044 
01045      newar = mri_read_dicom( new_fname );
01046 
01047    } else if( strstr(new_fname,".hdr") != NULL ||
01048               strstr(new_fname,".HDR") != NULL   ){  /* 05 Feb 2001 */
01049 
01050       newar = mri_read_analyze75( new_fname ) ;      /* ANALYZE .hdr/.img filepair */
01051 
01052    } else if( strstr(new_fname,".ima") != NULL ||
01053               strstr(new_fname,".IMA") != NULL   ){  /* 12 Mar 2001 */
01054 
01055       newar = mri_read_siemens( new_fname ) ;        /* Siemens file */
01056 
01057    } else if( strncmp(new_fname,"I.",2) == 0    ||  /* GE I.* files */
01058               strstr(new_fname,"/I.")   != NULL ||
01059               strstr(new_fname,".ppm")  != NULL ||  /* raw PPM or PGM files */
01060               strstr(new_fname,".pgm")  != NULL ||
01061               strstr(new_fname,".pnm")  != NULL ||
01062               strstr(new_fname,".PPM")  != NULL ||
01063               strstr(new_fname,".PNM")  != NULL ||
01064               strstr(new_fname,".PGM")  != NULL   ){ /* 05 Nov 2002 */
01065 
01066       newim = mri_read( new_fname ) ;      /* read from a 2D file with 1 slice */
01067 
01068       if ( newim == NULL )                 /* GEMS 4.x - 03 Jun 2003 [rickr] */
01069          newim = mri_read_ge4( new_fname ) ;
01070 
01071       if( newim != NULL ){
01072         INIT_IMARR(newar) ;
01073         ADDTO_IMARR(newar,newim) ;
01074       }
01075 
01076    } else if( strncmp(new_fname,"i.",2) == 0    ||  /* GEMS 4.x i.* files  */
01077               strstr(new_fname,"/i.")   != NULL ){  /* 03 Jun 2003 [rickr] */
01078 
01079       newim = mri_read_ge4( new_fname ) ;          /* 2D file with 1 slice */
01080 
01081       if( newim != NULL ){
01082         INIT_IMARR(newar) ;
01083         ADDTO_IMARR(newar,newim) ;
01084       }
01085 
01086    } else if( strstr(new_fname,".jpg" ) != NULL ||  /* various formats  */
01087               strstr(new_fname,".JPG" ) != NULL ||  /* that we convert  */
01088               strstr(new_fname,".jpeg") != NULL ||  /* to PPG/PGM using */
01089               strstr(new_fname,".JPEG") != NULL ||  /* external filters */
01090               strstr(new_fname,".gif" ) != NULL ||
01091               strstr(new_fname,".GIF" ) != NULL ||
01092               strstr(new_fname,".tif" ) != NULL ||
01093               strstr(new_fname,".TIF" ) != NULL ||
01094               strstr(new_fname,".tiff") != NULL ||
01095               strstr(new_fname,".TIFF") != NULL ||
01096               strstr(new_fname,".bmp" ) != NULL ||
01097               strstr(new_fname,".BMP" ) != NULL ||
01098               strstr(new_fname,".pbm" ) != NULL ||
01099               strstr(new_fname,".PBM" ) != NULL ||
01100               strstr(new_fname,".png" ) != NULL ||
01101               strstr(new_fname,".PNG" ) != NULL   ){ /* 22 Nov 2002 */
01102 
01103       newim = mri_read_stuff( new_fname ) ;
01104       if( newim != NULL ){
01105         INIT_IMARR(newar) ;
01106         ADDTO_IMARR(newar,newim) ;
01107       }
01108 
01109    } else if( strstr(new_fname,".mpg" ) != NULL ||  /* 03 Dec 2003 */
01110               strstr(new_fname,".MPG" ) != NULL ||  /* read MPEGs  */
01111               strstr(new_fname,".mpeg") != NULL ||
01112               strstr(new_fname,".MPEG") != NULL   ){
01113 
01114       newar = mri_read_mpeg( new_fname ) ;  /* cf. mri_read_mpeg.c */
01115    }
01116 
01117    /** failed to read anything?  try DICOM format (doesn't have a fixed suffix) **/
01118    /* 05 May 2003 added option to try DICOM last                    KRH          */
01119 
01120    if( newar == NULL ){
01121 
01122       if ( !AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
01123         newar = mri_read_dicom( new_fname ) ;  /* cf. mri_read_dicom.c */
01124       }
01125 
01126       /** if DICOM failed, try a 2D slice file, hope for the best **/
01127 
01128       if( newar == NULL ){
01129         newim = mri_read( new_fname ) ;
01130         if( newim == NULL ){ free(new_fname); RETURN( NULL ); }  /* give up */
01131         INIT_IMARR(newar) ;
01132         ADDTO_IMARR(newar,newim) ;
01133       }
01134 
01135       if ( (newar == NULL) && AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
01136         newar = mri_read_dicom( new_fname ) ;  /* cf. mri_read_dicom.c */
01137       }
01138    }
01139 
01140    free(new_fname) ;  /* done with the mangled filename */
01141 
01142    /* 07 Mar 2002: add fname to the images, if needed */
01143 
01144    if( newar != NULL && newar->num > 0 ){
01145      int ii ;
01146      for( ii=0 ; ii < newar->num ; ii++ ){
01147        newim = IMARR_SUBIM(newar,ii) ;
01148        if( newim != NULL && newim->fname == NULL )
01149           newim->fname = strdup(fname) ;
01150      }
01151    }
01152 
01153    RETURN( newar );
01154 }
01155 
01156 /*-----------------------------------------------------------------*/
01157 
01158 /*! Like mri_read_file(), but will only return 1 2D image.
01159 
01160     If the input file has more than 1 slice, or cannot be read,
01161     then NULL is returned.
01162 */
01163 
01164 MRI_IMAGE * mri_read_just_one( char * fname )
01165 {
01166    MRI_IMARR * imar ;
01167    MRI_IMAGE * im ;
01168    char * new_fname ;
01169 
01170 ENTRY("mri_read_just_one") ;
01171 
01172    new_fname = imsized_fname( fname ) ;
01173    if( new_fname == NULL ) RETURN( NULL );
01174 
01175    imar = mri_read_file( new_fname ) ; free(new_fname) ;
01176    if( imar == NULL ) RETURN( NULL );
01177    if( imar->num != 1 ){ DESTROY_IMARR(imar) ; RETURN( NULL ); }
01178    im = IMAGE_IN_IMARR(imar,0) ;
01179    FREE_IMARR(imar) ;
01180    RETURN( im );
01181 }
01182 
01183 /*-----------------------------------------------------------------
01184   return a count of how many 2D images will be read from this file
01185 -------------------------------------------------------------------*/
01186 
01187 /*! Return a count of how many 2D images are in a file.
01188 
01189     Used by to3d.c to figure out how many slices will be read
01190     later using mri_read_file().  Return value is 0 if the images
01191     can't be counted.  If you add a new file type to mri_read_file(),
01192     then you need to modify this function as well!
01193 */
01194 
01195 static int mri_imcount_analyze75( char * ) ;  /* prototype */
01196 static int mri_imcount_siemens( char * ) ;
01197 
01198 int mri_imcount( char * tname )
01199 {
01200    int hglobal , himage , nx , ny , nz , ngood ;
01201    char fname[256]="\0" ;
01202    char * new_fname ;
01203 
01204 ENTRY("mri_imcount") ;
01205 
01206    if( tname == NULL ) RETURN( 0 );
01207    new_fname = imsized_fname( tname ) ;
01208    if( new_fname == NULL ) RETURN( 0 );
01209 
01210    /*** a 3D filename ***/
01211 
01212    if( strlen(new_fname) > 9 && new_fname[0] == '3' && new_fname[1] == 'D' &&
01213        (new_fname[2] == ':' || new_fname[3] == ':') ){
01214                                /* check for ':', too   3 Jan 2005 [rickr] */
01215       switch( new_fname[2] ){
01216 
01217          default:
01218          case ':':
01219             ngood = sscanf( new_fname , "3D:%d:%d:%d:%d:%d:%s" ,
01220                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01221             break ;
01222 
01223          case 's':
01224             ngood = sscanf( new_fname , "3Ds:%d:%d:%d:%d:%d:%s" ,
01225                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01226             break ;
01227 
01228          case 'b':
01229             ngood = sscanf( new_fname , "3Db:%d:%d:%d:%d:%d:%s" ,
01230                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01231             break ;
01232 
01233          case 'f':
01234             ngood = sscanf( new_fname , "3Df:%d:%d:%d:%d:%d:%s" ,
01235                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01236             break ;
01237 
01238          case 'd':                                            /* 06 Feb 2003 */
01239             ngood = sscanf( new_fname , "3Dd:%d:%d:%d:%d:%d:%s" ,
01240                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01241             break ;
01242 
01243          case 'i':
01244             ngood = sscanf( new_fname , "3Di:%d:%d:%d:%d:%d:%s" ,
01245                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01246             break ;
01247 
01248          case 'c':
01249             ngood = sscanf( new_fname , "3Dc:%d:%d:%d:%d:%d:%s" ,
01250                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01251             break ;
01252 
01253          case 'r':
01254             ngood = sscanf( new_fname , "3Dr:%d:%d:%d:%d:%d:%s" ,
01255                             &hglobal , &himage , &nx , &ny , &nz , fname ) ;
01256             break ;
01257       }
01258 
01259       free( new_fname ) ;
01260       if( ngood < 6 || himage < 0 ||
01261           nx <= 0   || ny <= 0    || nz <= 0 ||
01262           strlen(fname) <= 0                       ) RETURN( 0 );
01263       else                                           RETURN( nz );
01264    }
01265 
01266    /*** a 3A filename ***/
01267 
01268    if( strlen(new_fname) > 9 &&
01269        new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
01270 
01271       switch( new_fname[2] ){
01272 
01273          default: ngood = 0 ; break ;
01274 
01275          case 's':
01276             ngood = sscanf( new_fname, "3As:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
01277             break ;
01278 
01279          case 'b':
01280             ngood = sscanf( new_fname, "3Ab:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
01281             break ;
01282 
01283          case 'f':
01284             ngood = sscanf( new_fname, "3Af:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
01285             break ;
01286       }
01287 
01288       free( new_fname ) ;
01289       if( ngood < 4 || nx <= 0 || ny <= 0 || nz <= 0 || strlen(fname) <= 0 ) RETURN( 0 );
01290       else                                                                   RETURN( nz );
01291    }
01292 
01293    /*** 05 Feb 2001: deal with ANALYZE .hdr files ***/
01294 
01295    if( strstr(new_fname,".hdr") != NULL ||
01296        strstr(new_fname,".HDR") != NULL   ){
01297 
01298       nz = mri_imcount_analyze75( new_fname ) ;
01299       if( nz > 0 ){ free(new_fname); RETURN(nz); }
01300    }
01301 
01302    if( strstr(new_fname,".ima") != NULL ||
01303        strstr(new_fname,".IMA") != NULL   ){        /* 12 Mar 2001 */
01304 
01305       nz = mri_imcount_siemens( new_fname ) ;
01306       if( nz > 0 ){ free(new_fname); RETURN(nz); }
01307    }
01308 
01309    if( strstr(new_fname,".mpg" ) != NULL ||  /* 03 Dec 2003 */
01310        strstr(new_fname,".MPG" ) != NULL ||
01311        strstr(new_fname,".mpeg") != NULL ||
01312        strstr(new_fname,".MPEG") != NULL   ){
01313 
01314       nz = mri_imcount_mpeg( new_fname ) ;
01315       if( nz > 0 ){ free(new_fname); RETURN(nz); }
01316    }
01317 
01318    /*** 19 Jul 2002: see if it is a DICOM file ***/
01319 
01320    mri_dicom_seterr(0) ;
01321    nz = mri_imcount_dicom( new_fname ) ;  /* cf. mri_read_dicom.c */
01322    mri_dicom_seterr(1) ;
01323    if( nz > 0 ){ free(new_fname); RETURN(nz); }
01324 
01325    /*** not recognized ***/
01326 
01327    free(new_fname) ; RETURN(1) ;    /* assume it has 1 image in it, somewhere */
01328 }
01329 
01330 /*--------------------------------------------------------------*/
01331 
01332 /*! Like mri_read_file(), but returns images from many files.
01333 
01334     \param nf = Number of file names
01335     \param fn = Array of file name strings
01336     \return An array of 2D images (NULL if nothing was found)
01337 
01338     Added 07 Mar 1995
01339 */
01340 
01341 MRI_IMARR * mri_read_many_files( int nf , char * fn[] )
01342 {
01343    MRI_IMARR * newar , * outar ;
01344    int kf , ii ;
01345 
01346 ENTRY("mri_read_many_files") ;
01347 
01348    if( nf <= 0 ) RETURN( NULL );  /* no inputs! */
01349    INIT_IMARR(outar) ;          /* initialize output array */
01350 
01351    for( kf=0 ; kf < nf ; kf++ ){
01352       newar = mri_read_file( fn[kf] ) ;  /* read all images in this file */
01353 
01354       if( newar == NULL ){  /* none?  flush the output array! */
01355          fprintf(stderr,"cannot read images from file %s\n",fn[kf]) ;
01356          for( ii=0 ; ii < outar->num ; ii++ ) mri_free(outar->imarr[ii]) ;
01357          FREE_IMARR(outar) ;
01358          RETURN( NULL );
01359       }
01360 
01361       for( ii=0 ; ii < newar->num ; ii++ )  /* move images to output array */
01362          ADDTO_IMARR( outar , newar->imarr[ii] ) ;
01363 
01364       FREE_IMARR(newar) ;  /* don't need this no more */
01365    }
01366    RETURN( outar );
01367 }
01368 
01369 /*---------------------------------------------------------------*/
01370 
01371 /*! Read a raw PPM file into 3 byte-valued MRI_IMAGEs.
01372 
01373     \date 16 May 1995
01374 */
01375 
01376 MRI_IMARR * mri_read_ppm3( char * fname )
01377 {
01378    int ch , nch , nx,ny,maxval , length , npix,ii ;
01379    char buf[512] ;
01380    MRI_IMAGE *rim , *gim , *bim ;
01381    MRI_IMARR * outar ;
01382    FILE * imfile ;
01383    byte * rby , * gby , * bby , * rgby ;
01384 
01385 ENTRY("mri_read_ppm3") ;
01386 
01387    /*** open input file ***/
01388 
01389    imfile = fopen( fname , "r" ) ;
01390    if( imfile == NULL ){
01391       fprintf(stderr,"couldn't open file %s in mri_read_ppm3\n",fname); RETURN(NULL) ;
01392    }
01393 
01394    /*** check if a raw PPM file ***/
01395 
01396    ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; RETURN(NULL); }
01397    ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; RETURN(NULL); }
01398 
01399    /* magic P6 found, so read numbers in header */
01400 
01401    ch = getc(imfile) ;
01402 
01403    NUMSCAN(nx)     ; if( nx     <= 0 )   { fclose(imfile) ; RETURN(NULL); }
01404    NUMSCAN(ny)     ; if( ny     <= 0 )   { fclose(imfile) ; RETURN(NULL); }
01405    NUMSCAN(maxval) ; if( maxval <= 0 ||
01406                          maxval >  255 ) { fclose(imfile) ; RETURN(NULL); }
01407 
01408    /*** create output images and workspace array ***/
01409 
01410    rim = mri_new( nx , ny , MRI_byte ) ; rby = mri_data_pointer( rim ) ;
01411    gim = mri_new( nx , ny , MRI_byte ) ; gby = mri_data_pointer( gim ) ;
01412    bim = mri_new( nx , ny , MRI_byte ) ; bby = mri_data_pointer( bim ) ;
01413 
01414    sprintf(buf,"%s#R",fname) ; mri_add_name( buf , rim ) ;
01415    sprintf(buf,"%s#G",fname) ; mri_add_name( buf , gim ) ;
01416    sprintf(buf,"%s#B",fname) ; mri_add_name( buf , bim ) ;
01417 
01418    rgby = (byte *) malloc( sizeof(byte) * 3*nx*ny ) ;
01419    if( rgby == NULL ){
01420       fprintf(stderr,"couldn't malloc workspace in mri_read_ppm3!\n") ; EXIT(1) ;
01421    }
01422 
01423    /*** read all data into workspace array ***/
01424 
01425    length = fread( rgby , sizeof(byte) , 3*nx*ny , imfile ) ;
01426    fclose( imfile ) ;
01427 
01428    if( length != 3*nx*ny ){
01429       free(rgby) ; mri_free(rim) ; mri_free(gim) ; mri_free(bim) ;
01430       fprintf(stderr,"couldn't read data from file %s in mri_read_ppm3\n",fname) ;
01431       RETURN(NULL);
01432    }
01433 
01434    /*** put data from workspace array into output images ***/
01435 
01436    npix = nx*ny ;
01437    for( ii=0 ; ii < npix ; ii++ ){
01438       rby[ii] = rgby[3*ii  ] ;
01439       gby[ii] = rgby[3*ii+1] ;
01440       bby[ii] = rgby[3*ii+2] ;
01441    }
01442    free( rgby ) ;
01443 
01444    /*** create output image array ***/
01445 
01446    INIT_IMARR(outar) ;
01447    ADDTO_IMARR( outar , rim ) ;
01448    ADDTO_IMARR( outar , gim ) ;
01449    ADDTO_IMARR( outar , bim ) ;
01450    RETURN(outar);
01451 }
01452 
01453 /*-----------------------------------------------------------------
01454    routines added 1 Oct 1995
01455 -------------------------------------------------------------------*/
01456 
01457 /*! Read 1 2D image, then "nsize" it - make it a power of 2 in sizes.
01458 
01459     This was developed in the days when FD/FD2/fim ruled the world, and
01460     those programs (AJ's legacy) only deal with square images that are
01461     a power of 2 in size.
01462     \date 01 Oct 1995
01463 */
01464 
01465 MRI_IMAGE *mri_read_nsize( char * fname )
01466 {
01467    MRI_IMARR *imar ;
01468    MRI_IMAGE *imout ;
01469 
01470    imar = mri_read_file( fname ) ;
01471    if( imar == NULL ) return NULL ;
01472    if( imar->num != 1 ){ DESTROY_IMARR(imar) ; return NULL ; }
01473 
01474    imout = mri_nsize( IMAGE_IN_IMARR(imar,0) ) ;
01475    mri_add_name( IMAGE_IN_IMARR(imar,0)->name , imout ) ;
01476 
01477    DESTROY_IMARR(imar) ;
01478    return imout ;
01479 }
01480 
01481 /*! Read many 2D images from many files. */
01482 
01483 MRI_IMARR *mri_read_many_nsize( int nf , char * fn[] )
01484 {
01485    MRI_IMARR * newar , * outar ;
01486    MRI_IMAGE * im ;
01487    int ii ;
01488 
01489    newar = mri_read_many_files( nf , fn ) ;
01490    if( newar == NULL ) return NULL ;
01491 
01492    INIT_IMARR(outar) ;
01493    for( ii=0 ; ii < newar->num ; ii++ ){
01494       im = mri_nsize( IMAGE_IN_IMARR(newar,ii) ) ;
01495       mri_add_name( IMAGE_IN_IMARR(newar,ii)->name , im ) ;
01496       ADDTO_IMARR(outar,im) ;
01497       mri_free( IMAGE_IN_IMARR(newar,ii) ) ;
01498    }
01499    FREE_IMARR(newar) ;
01500    return outar ;
01501 }
01502 
01503 /*------------------------------------------------------------------------*/
01504 
01505 /*! Set up MCW_SIZE_# database for input.
01506 
01507     This implements the facility for the user to define MCW_IMSIZE_1
01508     (or AFNI_IMSIZE_1) et cetera, for pre-defining a relationship between
01509     a file size in bytes and a 3D: prefix.  This function is only called
01510     once to setup the table.
01511     \date 07 Nov 95
01512 */
01513 
01514 void init_MCW_sizes(void)
01515 {
01516    int num , count ;
01517    char ename[32] ;
01518    char * str ;
01519 
01520    if( MCW_imsize_good >= 0 ) return ;
01521 
01522    MCW_imsize_good = 0 ;
01523 
01524    for( num=0 ; num < MAX_MCW_IMSIZE ; num++ ){ /* look for environment string */
01525 
01526       imsize[num].size = -1 ;
01527 
01528       /* try to find environment variable with the num-th name */
01529 
01530       sprintf( ename , "AFNI_IMSIZE_%d" , num+1 ) ;
01531       str = my_getenv( ename ) ;
01532 
01533       if( str == NULL ){
01534          sprintf( ename , "MCW_IMSIZE_%d" , num+1 ) ;
01535          str = my_getenv( ename ) ;
01536          if( str == NULL ) continue ;
01537       }
01538 
01539       imsize[num].prefix = (char *) malloc( sizeof(char) * strlen(str) ) ;
01540       if( imsize[num].prefix == NULL ){
01541          fprintf(stderr,"\n*** Can't malloc in init_MCW_sizes! ***\a\n");
01542          EXIT(1) ;
01543       }
01544 
01545       if( str[0] != '%' ){  /* e.g., 16096=3D:-1:0:64:64:1: */
01546 
01547          imsize[num].head = -1 ;
01548          count = sscanf( str , "%d=%s" , &(imsize[num].size) , imsize[num].prefix ) ;
01549          if( count != 2 || imsize[num].size < 2 || strlen(imsize[num].prefix) < 2 ){
01550             free( imsize[num].prefix ) ;
01551             fprintf(stderr,"bad environment %s = %s\n" ,
01552                     ename , str ) ;
01553          }
01554 
01555       } else {              /* e.g., %16096+0=3D:0:7904:64:64: */
01556 
01557          count = sscanf( str+1 , "%d+%d=%s" ,
01558                          &(imsize[num].size) , &(imsize[num].head) , imsize[num].prefix ) ;
01559 
01560          if( count != 3 || imsize[num].size < 2 ||
01561              imsize[num].head < 0 || strlen(imsize[num].prefix) < 2 ){
01562 
01563             free( imsize[num].prefix ) ;
01564             fprintf(stderr,"bad environment %s = %s\n" ,
01565                     ename , str ) ;
01566          }
01567       }
01568 
01569       MCW_imsize_good ++ ;
01570    }
01571 
01572    return ;
01573 }
01574 
01575 /*------------------------------------------------------------------------------*/
01576 /*! My version of strdup(), which won't fail if the input is NULL. */
01577 
01578 char * my_strdup( char * str )
01579 {
01580    char * new_str ;
01581    if( str == NULL ) return NULL ;
01582    new_str = (char *) malloc( sizeof(char) * (strlen(str)+1) ) ;
01583    if( new_str != NULL ) strcpy( new_str , str ) ;
01584    return new_str ;
01585 }
01586 
01587 /*------------------------------------------------------------------------------*/
01588 
01589 /*! Check if a filesize fits an MCW_IMSIZE setup.
01590 
01591     \param fname = Filename
01592     \return A new "filename" with 3D header attached if it fits.
01593             If not, return a copy of the filename.  In any case the
01594             returned string should be free()-d when it is no longer needed.
01595 */
01596 
01597 char * imsized_fname( char * fname )
01598 {
01599    int num , lll ;
01600    long len ;
01601    char * new_name ;
01602 
01603    init_MCW_sizes() ;
01604    if( MCW_imsize_good == 0 ){
01605       new_name = my_strdup(fname) ;  /* nothing to fit */
01606       return new_name ;              /* --> return copy of old name */
01607    }
01608 
01609    len = mri_filesize( fname ) ;
01610    if( len <= 0 ){
01611       new_name = my_strdup(fname) ;  /* not an existing filename */
01612       return new_name ;              /* --> return copy of old name */
01613    }
01614 
01615    for( num=0 ; num < MAX_MCW_IMSIZE ; num++ ){     /* check each possibility */
01616 
01617       if( imsize[num].size <= 0 ) continue ;        /* skip to next one */
01618 
01619       if( imsize[num].head < 0 && len == imsize[num].size ){  /* fixed size fit */
01620 
01621          lll = strlen(fname) + strlen(imsize[num].prefix) + 4 ;
01622          new_name = (char *) malloc( sizeof(char) * lll ) ;
01623          if( new_name == NULL ){
01624             fprintf(stderr,"\n*** Can't malloc in imsized_fname! ***\a\n");
01625             EXIT(1) ;
01626          }
01627          sprintf( new_name , "%s%s" , imsize[num].prefix , fname ) ;
01628          return new_name ;
01629 
01630       } else if( (len-imsize[num].head) % imsize[num].size == 0 ){
01631          int count = (len-imsize[num].head) / imsize[num].size ;
01632 
01633          if( count < 1 ) continue ;  /* skip to next one */
01634 
01635          lll = strlen(fname) + strlen(imsize[num].prefix) + 32 ;
01636          new_name = (char *) malloc( sizeof(char) * lll ) ;
01637          if( new_name == NULL ){
01638             fprintf(stderr,"\n*** Can't malloc in imsized_fname! ***\a\n");
01639             EXIT(1) ;
01640          }
01641          sprintf( new_name , "%s%d:%s" , imsize[num].prefix , count , fname ) ;
01642          return new_name ;
01643       }
01644 
01645    }
01646 
01647    new_name = my_strdup(fname) ;  /* no fit --> return copy of old name */
01648    return new_name ;
01649 }
01650 
01651 /*------------------------------------------------------------------------*/
01652 /*! Return the size of a file in bytes.
01653 
01654   \param pathname = input filename
01655   \return File length if file exists; -1 if it doesn't.
01656   \see THD_filesize() in thd_filestuff.c.
01657 */
01658 
01659 long mri_filesize( char * pathname )
01660 {
01661    static struct stat buf ;
01662    int ii ;
01663 
01664    if( pathname == NULL ) return -1 ;
01665    ii = stat( pathname , &buf ) ; if( ii != 0 ) return -1 ;
01666    return buf.st_size ;
01667 }
01668 
01669 /*---------------------------------------------------------------*/
01670 
01671 /*! Read the header from PPM file and return its info.
01672 
01673   \param fname = file name
01674   \return *nx and *ny are set to the image dimensions;
01675           if they are set to 0, something bad happened
01676           (e.g., the file isn't a PPM file, or doesn't exist).
01677   \date 17 Sep 2001
01678 */
01679 
01680 void mri_read_ppm_header( char *fname , int *nx, int *ny )
01681 {
01682    FILE *imfile ;
01683    int ch , nch , nxx,nyy ;
01684    char buf[256] ;
01685 
01686 ENTRY("mri_read_ppm_header") ;
01687 
01688    if( fname == NULL || nx == NULL || ny == NULL ) EXRETURN ;
01689 
01690    *nx = *ny = 0 ;  /* default returns */
01691 
01692    /*** open input file ***/
01693 
01694    imfile = fopen( fname , "r" ) ; if( imfile == NULL ) EXRETURN ;
01695 
01696    /*** check if a raw PPM file ***/
01697 
01698    ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; EXRETURN ; }
01699    ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; EXRETURN ; }
01700 
01701    /* magic P6 found, so read numbers in header */
01702 
01703    ch = getc(imfile) ;
01704 
01705    NUMSCAN(nxx) ; if( nxx <= 0 ){ fclose(imfile) ; EXRETURN ; }
01706    NUMSCAN(nyy) ; if( nyy <= 0 ){ fclose(imfile) ; EXRETURN ; }
01707 
01708    /* return dimensions */
01709 
01710    fclose(imfile) ; *nx = nxx ; *ny = nyy ; EXRETURN ;
01711 }
01712 
01713 /*---------------------------------------------------------------*/
01714 
01715 /*! Reads a raw PPM file into 1 2D MRI_rgb-valued image.
01716 
01717    \param fname = Image filename
01718    \return An MRI_IMAGE if things worked OK; NULL if not
01719    \date 13 May 1996
01720 */
01721 
01722 MRI_IMAGE * mri_read_ppm( char * fname )
01723 {
01724    int ch , nch , nx,ny,maxval , length ;
01725    MRI_IMAGE * rgbim ;
01726    FILE      * imfile ;
01727    byte      * rgby ;
01728    char        buf[256] ;
01729 
01730 ENTRY("mri_read_ppm") ;
01731 
01732    /*** open input file ***/
01733 
01734    imfile = fopen( fname , "r" ) ;
01735    if( imfile == NULL ) RETURN(NULL);
01736 
01737    /*** check if a raw PPM file ***/
01738 
01739    ch = getc( imfile ) ; if( ch != 'P' ) { fclose(imfile) ; RETURN(NULL); }
01740    ch = getc( imfile ) ; if( ch != '6' ) { fclose(imfile) ; RETURN(NULL); }
01741 
01742    /* magic P6 found, so read numbers in header */
01743 
01744    ch = getc(imfile) ;
01745 
01746    NUMSCAN(nx)    ; if( nx     <= 0 )  { fclose(imfile); RETURN(NULL); }
01747    NUMSCAN(ny)    ; if( ny     <= 0 )  { fclose(imfile); RETURN(NULL); }
01748    NUMSCAN(maxval); if( maxval <= 0 ||
01749                         maxval >  255 ){ fclose(imfile); RETURN(NULL); }
01750 
01751    /*** create output image ***/
01752 
01753    rgbim = mri_new( nx , ny , MRI_rgb ) ; mri_add_name( fname , rgbim ) ;
01754    rgby  = MRI_RGB_PTR(rgbim) ;
01755 
01756    /*** read all data into image array */
01757 
01758    length = fread( rgby , sizeof(byte) , 3*nx*ny , imfile ) ;
01759    fclose( imfile ) ;
01760 
01761    if( length != 3*nx*ny ){ mri_free(rgbim) ; RETURN(NULL) ; }
01762 
01763    /* 17 Sep 2001: scale to maxval=255, if needed */
01764 
01765    if( maxval < 255 ){
01766       int ii ; float fac = 255.4/maxval ;
01767       for( ii=0 ; ii < 3*nx*ny ; ii++ ) rgby[ii] = (byte)( rgby[ii]*fac ) ;
01768    }
01769 
01770    RETURN(rgbim) ;
01771 }
01772 /*---------------------------------------------------------------*/
01773 
01774 /*! Length of line buffer for mri_read_ascii() */
01775 /* rcr - improve this */
01776 #define LBUF 2524288  /* 08 Jul 2004: increased to 512K from 64K */
01777 
01778 /*! Free a buffer and set it to NULL */
01779 #define FRB(b) do{ if( (b)!=NULL ){free((b)); (b)=NULL;} }while(0)
01780 
01781 #undef USE_LASTBUF
01782 
01783 /*---------------------------------------------------------------*/
01784 /*! [20 Jun 2002] Like fgets, but also
01785      - skips blank or comment lines
01786      - skips leading and trailing whitespace
01787      - catenates lines that end in '\' (replacing '\' with ' ')
01788      - returns duplicate of last line if first 2
01789         nonblank input characters are "" [20 Jul 2004]
01790 -----------------------------------------------------------------*/
01791 
01792 static char * my_fgets( char *buf , int size , FILE *fts )
01793 {
01794    char *ptr ;
01795    int nbuf , ll,ii , cflag ;
01796    static char *qbuf=NULL ;
01797 
01798 #ifdef USE_LASTBUF
01799    static char *lastbuf = NULL ;   /* 20 Jul 2004 */
01800    static int  nlastbuf = 0 ;
01801 
01802    if( buf == NULL && lastbuf != NULL ){    /* 20 Jul 2004 */
01803      free((void *)lastbuf); lastbuf = NULL; nlastbuf = 0 ;
01804    }
01805 #endif
01806 
01807    if( buf == NULL && qbuf != NULL ){ free((void *)qbuf); qbuf = NULL; }
01808 
01809    if( buf == NULL || size < 1 || fts == NULL ) return NULL ;
01810 
01811    if( qbuf == NULL ) qbuf = AFMALL(char, LBUF) ;  /* 1st time in */
01812 
01813    nbuf  = 0 ;  /* num bytes stored in buf so far */
01814    cflag = 0 ;  /* flag if we're catenating lines */
01815 
01816    while(1){   /* loop and read lines, creating a logical line */
01817 
01818      ptr = fgets( qbuf , LBUF , fts ) ; /* read next whole line */
01819 
01820      if( ptr == NULL ) break ;          /* must be end-of-file */
01821 
01822      /* skip leading whitespace */
01823 
01824      for( ; *ptr != '\0' && isspace(*ptr) ; ptr++ ) ; /* nada */
01825 
01826      /* skip entirely blank lines, unless we are catenating */
01827 
01828      if( *ptr == '\0' ){ if(cflag) break; else continue; }
01829 
01830 #ifdef USE_LASTBUF
01831      /* if a duplicate is requested, return it now [20 Jul 2004] */
01832 
01833      if( *ptr == '"' && *(ptr+1) == '"' && nlastbuf > 0 && nbuf == 0 ){
01834        ll = strlen(lastbuf) ; if( ll >= size ) ll = size-1 ;
01835        memcpy(buf,lastbuf,ll-1) ; buf[ll] = '\0' ;
01836        return buf ;
01837      }
01838 #endif
01839 
01840      /* skip comment lines (even if we are catenating) */
01841 
01842      if( *ptr == '#' || (*ptr == '/' && *(ptr+1) == '/') ) continue ;
01843 
01844      /* strip trailing whitespace */
01845 
01846      ll = strlen(ptr) ;                                  /* will be > 0 */
01847      for( ii=ll-1 ; isspace(ptr[ii]) && ii > 0 ; ii-- )  /* blank => NUL */
01848        ptr[ii] = '\0' ;
01849 
01850      ll = strlen(ptr) ;                 /* number of chars left */
01851      if( ll == 0 ) continue ;           /* should not happen */
01852 
01853      cflag = (ptr[ll-1] == '\\') ;      /* catenate next line? */
01854      if( cflag ) ptr[ll-1] = ' ' ;      /* replace '\' with ' ' */
01855 
01856      /* now copy what's left (ll+1 bytes) at tail of output buffer */
01857 
01858      if( nbuf+ll+1 > size ){   /* too much for output buffer? */
01859        ll = size - (nbuf+1) ;
01860        if( ll <= 0 ) break ;   /* should not happen */
01861      }
01862 
01863      memcpy(buf+nbuf,ptr,ll+1) ; nbuf += ll ;
01864      if( !cflag ) break ;
01865 
01866    } /* loop to get next line if catenation is turned on */
01867 
01868 #ifdef LASTBUF
01869    /* make a copy of result in lastbuf [20 Jul 2004] */
01870 
01871    ll = strlen(buf) ;
01872    if( ll+1 > nlastbuf ){
01873      nlastbuf = ll+2 ; lastbuf = (char *)realloc((void *)lastbuf,nlastbuf) ;
01874    }
01875    memcpy(lastbuf,buf,ll+1) ;
01876 #endif
01877 
01878    /* and we is done */
01879 
01880    if( nbuf > 0 ) return buf ;      /* return what we read already */
01881    return NULL ;                    /* signal of failure get data  */
01882 }
01883 
01884 /*--------------------------------------------------------------*/
01885 static float lbfill = 0.0 ;  /* 10 Aug 2004 */
01886 
01887 /*--------------------------------------------------------------*/
01888 /*! Decode a line buffer into an array of floats.               */
01889 
01890 static floatvec * decode_linebuf( char *buf )  /* 20 Jul 2004 */
01891 {
01892    floatvec *fv=NULL ;
01893    int blen, bpos, ncol, ii, count ;
01894    char sep, vbuf[64] , *cpt ;
01895    float val ;
01896 
01897    if( buf == NULL || *buf == '\0' ) return fv ;
01898 
01899    blen = strlen(buf) ;
01900    ncol = 0 ;
01901 
01902    /* convert commas to blanks */
01903 
01904    for( ii=0 ; ii < blen ; ii++ ) if( buf[ii] == ',' ) buf[ii] = ' ' ;
01905 
01906    fv = (floatvec *)malloc(sizeof(floatvec)) ;
01907    fv->nar = 0 ;
01908    fv->ar  = (float *)NULL ;
01909 
01910    for( bpos=0 ; bpos < blen ; ){
01911      /* skip to next nonblank character */
01912 
01913      for( ; bpos < blen && (isspace(buf[bpos])||buf[bpos]==',') ; bpos++ ) ; /* nada */
01914      if( bpos == blen ) break ;    /* end of line */
01915 
01916      sscanf( buf+bpos , "%63s" , vbuf ) ;
01917 
01918      val = 0.0 ; count = 1 ;
01919      if( vbuf[0] == '*' ){    /* 10 Aug 2004 */
01920        val = lbfill ;
01921      } else if( (cpt=strchr(vbuf,'@')) != NULL ){
01922        sscanf( vbuf , "%d%c%f" , &count , &sep , &val ) ;
01923        if( count < 1 ) count = 1 ;
01924        if( *(cpt+1) == '*' ) val = lbfill ;  /* 10 Aug 2004 */
01925      } else {
01926        sscanf( vbuf , "%f" , &val ) ;
01927      }
01928 
01929      fv->ar = (float *)realloc( (void *)fv->ar , sizeof(float)*(fv->nar+count) ) ;
01930      for( ii=0 ; ii < count ; ii++ ) fv->ar[ii+fv->nar] = val ;
01931      fv->nar += count ;
01932      bpos += strlen(vbuf) ;
01933    }
01934 
01935    if( fv->nar == 0 ){ KILL_floatvec(fv); fv = NULL; }
01936    return fv ;
01937 }
01938 
01939 /*---------------------------------------------------------------*/
01940 
01941 /*! Increment for time series array size for mri_read_ascii() */
01942 #define INC_TSARSIZE 128
01943 
01944 /*! Read an array of ASCII numbers into a 1D or 2D image.
01945 
01946   \param fname = input filename
01947   \return Pointer to MRI_IMAGE (in MRI_float) format if things
01948           are cool; NULL if not.
01949   \date Jun 1996
01950 
01951   Example input:
01952      - Line 1:  3 4 6
01953      - Line 2:  2 2 2
01954      - Line 3:  7 2 1
01955      - Line 4:  9 9 6
01956   This produces an image with nx=3 and ny=4.  The first row
01957   is read to determine nx; all subsequent rows must have nx
01958   values.  A line whose very first character is a '#' will
01959   be skipped as a comment.  A line with no characters (just
01960   the '\n') will also be skipped.
01961 
01962   20 Jun 2002: modified to use my_fgets() instead of fgets().
01963 */
01964 
01965 MRI_IMAGE * mri_read_ascii( char * fname )
01966 {
01967    MRI_IMAGE * outim ;
01968    int ii,jj,val , used_tsar , alloc_tsar ;
01969    float * tsar ;
01970    float ftemp ;
01971    FILE * fts ;
01972    char * ptr ;
01973    int  ncol , bpos , blen , nrow ;
01974    static char *buf=NULL ;            /* 20 Jun 2002: make a ptr */
01975 
01976    floatvec *fvec ;                   /* 20 Jul 2004 */
01977    int incts ;
01978 
01979 ENTRY("mri_read_ascii") ;
01980 
01981    if( fname == NULL || fname[0] == '\0' ) RETURN(NULL) ;
01982 
01983    if( strncmp(fname,"1D:",3) == 0 ){         /* 28 Apr 2003 */
01984      MRI_IMAGE *qim = mri_1D_fromstring( fname+3 ) ;
01985      if( qim != NULL ){
01986        outim = mri_transpose(qim); mri_free(qim); RETURN(outim);
01987      }
01988    }
01989 
01990    fts = fopen( fname , "r" ); if( fts == NULL ) RETURN(NULL);
01991 
01992    if( buf == NULL ) buf = AFMALL(char, LBUF) ; /* create buffer */
01993 
01994    /** step 1: read in the first line and see how many numbers are in it
01995                (skipping lines that are comments or entirely blank)     */
01996 
01997    (void) my_fgets( NULL , 0 , NULL ) ;  /* reset [20 Jul 2004] */
01998    ptr = my_fgets( buf , LBUF , fts ) ;
01999    if( ptr==NULL || *ptr=='\0' ){ FRB(buf); fclose(fts); RETURN(NULL); }  /* bad read? */
02000 
02001    lbfill = 0.0f ;                          /* 10 Aug 2004 */
02002 
02003    fvec = decode_linebuf( buf ) ;           /* 20 Jul 2004 */
02004    if( fvec == NULL || fvec->nar == 0 ){
02005      if( fvec != NULL ) KILL_floatvec(fvec) ;
02006      FRB(buf); fclose(fts); RETURN(NULL);
02007    }
02008    ncol = fvec->nar ; KILL_floatvec(fvec) ;
02009 
02010    /** At this point, ncol is the number of floats to be read from each line **/
02011 
02012    rewind( fts ) ;  /* will start over */
02013 
02014    incts      = MAX(INC_TSARSIZE,ncol) ;
02015    used_tsar  = 0 ;
02016    alloc_tsar = incts ;
02017    tsar       = (float *) malloc( sizeof(float) * alloc_tsar ) ;
02018    if( tsar == NULL ){
02019       fprintf(stderr,"\n*** malloc error in mri_read_ascii ***\n"); EXIT(1);
02020    }
02021 
02022    /** read lines, convert to floats, store **/
02023 
02024    nrow = 0 ;
02025    while( 1 ){
02026      ptr = my_fgets( buf , LBUF , fts ) ;  /* read */
02027      if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
02028 
02029      fvec = decode_linebuf( buf ) ;
02030      if( fvec == NULL ) break ;
02031      if( fvec->nar == 0 ){ KILL_floatvec(fvec); break; }
02032 
02033      if( used_tsar + ncol >= alloc_tsar ){
02034         alloc_tsar += incts ;
02035         tsar        = (float *)realloc( (void *)tsar,sizeof(float)*alloc_tsar );
02036         if( tsar == NULL ){
02037           fprintf(stderr,"\n*** realloc error in mri_read_ascii ***\n"); EXIT(1);
02038         }
02039      }
02040      for( ii=0 ; ii < fvec->nar && ii < ncol ; ii++ )
02041        tsar[used_tsar+ii] = fvec->ar[ii] ;
02042      for( ; ii < ncol ; ii++ )
02043        tsar[used_tsar+ii] = 0.0 ;
02044      used_tsar += ncol ;
02045      KILL_floatvec(fvec) ;
02046 
02047      nrow++ ;                  /* got one more complete row! */
02048    }
02049    fclose( fts ) ; /* finished with this file! */
02050    (void) my_fgets( NULL , 0 , NULL ) ;  /* reset [20 Jul 2004] */
02051 
02052    if( used_tsar <= 1 ){ FRB(buf); free(tsar); RETURN(NULL); }
02053 
02054    tsar = (float *) realloc( tsar , sizeof(float) * used_tsar ) ;
02055    if( tsar == NULL ){
02056       fprintf(stderr,"\n*** final realloc error in mri_read_ascii ***\n"); EXIT(1);
02057    }
02058 
02059    outim = mri_new_vol_empty( ncol , nrow , 1 , MRI_float ) ;
02060    mri_fix_data_pointer( tsar , outim ) ;
02061    mri_add_name( fname , outim ) ;
02062 
02063    FRB(buf) ; RETURN(outim) ;
02064 }
02065 
02066 /*---------------------------------------------------------------------------*/
02067 
02068 /*! Read an ASCII file as columns, transpose to rows, allow column selectors.
02069 
02070   \param fname = Input filename (max of 255 characters)
02071   \return Pointer to MRI_IMAGE if all went well; NULL if not.
02072   \date 16 Nov 1999
02073 
02074   This function builds on mri_read_ascii() in two ways:
02075     - the input is transposed to rows (so that a 1x100 file becomes a 100x1 image)
02076     - column selectors are allowed in fname
02077 */
02078 
02079 MRI_IMAGE * mri_read_1D( char * fname )
02080 {
02081    MRI_IMAGE *inim , *outim , *flim ;
02082    char dname[512] , *cpt , *dpt ;
02083    int ii,jj,nx,ny,nts , *ivlist , *ivl , *sslist ;
02084    float *far , *oar ;
02085 
02086 ENTRY("mri_read_1D") ;
02087 
02088    if( fname == NULL || fname[0] == '\0' || strlen(fname) > 511 ) RETURN(NULL) ;
02089 
02090    if( strncmp(fname,"1D:",3) == 0 ){       /* 28 Apr 2003 */
02091      return mri_1D_fromstring( fname+3 ) ;
02092    }
02093 
02094    /*-- split filename and subvector list --*/
02095 
02096    cpt = strstr(fname,"[") ;
02097    dpt = strstr(fname,"{") ;            /* 30 Apr 2003: subsampling list */
02098 
02099    if( cpt == fname || dpt == fname ){  /* can't be at start of filename! */
02100       fprintf(stderr,"*** Illegal filename in mri_read_1D: %s\n",fname) ;
02101       RETURN(NULL) ;
02102    } else {                             /* got a subvector list */
02103       strcpy( dname , fname ) ;
02104       if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; }
02105       if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; }
02106    }
02107 
02108    /*-- read file in --*/
02109 
02110    inim = mri_read_ascii(dname) ;
02111    if( inim == NULL ) RETURN(NULL) ;
02112    flim = mri_transpose(inim) ; mri_free(inim) ;
02113 
02114    /*-- get the subvector and subsampling lists, if any --*/
02115 
02116    nx = flim->nx ; ny = flim->ny ;
02117 
02118    ivlist = MCW_get_intlist( ny , cpt ) ;   /* subvector list */
02119    sslist = MCW_get_intlist( nx , dpt ) ;   /* subsampling list */
02120 
02121    /* if have subvector list, extract those rows into a new image */
02122 
02123    if( ivlist != NULL && ivlist[0] > 0 ){
02124      nts = ivlist[0] ;                         /* number of subvectors */
02125      ivl = ivlist + 1 ;                        /* start of array of subvectors */
02126 
02127      for( ii=0 ; ii < nts ; ii++ ){            /* check them out */
02128        if( ivl[ii] < 0 || ivl[ii] >= ny ){
02129          fprintf(stderr,"*** Out-of-range subvector [list] in mri_read_1D: %s\n",fname) ;
02130          mri_free(flim) ; free(ivlist) ; RETURN(NULL) ;
02131        }
02132      }
02133 
02134      outim = mri_new( nx , nts , MRI_float ) ; /* make output image */
02135      far   = MRI_FLOAT_PTR( flim ) ;
02136      oar   = MRI_FLOAT_PTR( outim ) ;
02137 
02138      for( ii=0 ; ii < nts ; ii++ )             /* copy desired rows */
02139        memcpy( oar + ii*nx , far + ivl[ii]*nx , sizeof(float)*nx ) ;
02140 
02141      mri_free(flim); free(ivlist); flim = outim; ny = nts;
02142    }
02143 
02144    /* if have subsampling list, extract those columns into a new image */
02145 
02146    if( sslist != NULL && sslist[0] > 0 ){
02147      nts = sslist[0] ;                         /* number of columns to get */
02148      ivl = sslist + 1 ;                        /* start of array of column indexes */
02149 
02150      for( ii=0 ; ii < nts ; ii++ ){            /* check them out */
02151        if( ivl[ii] < 0 || ivl[ii] >= nx ){
02152          fprintf(stderr,"*** Out-of-range subsampling {list} in mri_read_1D: %s\n",fname) ;
02153          mri_free(flim) ; free(sslist) ; RETURN(NULL) ;
02154        }
02155      }
02156 
02157      outim = mri_new( nts , ny , MRI_float ) ; /* make output image */
02158      far   = MRI_FLOAT_PTR( flim ) ;
02159      oar   = MRI_FLOAT_PTR( outim ) ;
02160 
02161      for( ii=0 ; ii < nts ; ii++ )             /* copy desired columns */
02162        for( jj=0 ; jj < ny ; jj++ )
02163          oar[ii+jj*nts] = far[ivl[ii]+jj*nx] ;
02164 
02165      mri_free(flim); free(sslist); flim = outim;
02166    }
02167 
02168    RETURN(flim) ;
02169 }
02170 
02171 /*-----------------------------------------------------------------------------------*/
02172 
02173 MRI_IMAGE * mri_read_ascii_ragged( char *fname , float filler )  /* 28 Jul 2004 - ragged rows */
02174 {
02175    MRI_IMAGE *outim ;
02176    int ii,jj , ncol,nrow ;
02177    float *tsar ;
02178    FILE *fts ;
02179    char *ptr ;
02180    static char *buf=NULL ;
02181    floatvec *fvec ;
02182 
02183 ENTRY("mri_read_ascii_ragged") ;
02184 
02185    if( fname == NULL || *fname == '\0' ){ FRB(buf); RETURN(NULL); }
02186 
02187    fts = fopen( fname , "r" ); if( fts == NULL ){ FRB(buf); RETURN(NULL); }
02188 
02189    if( buf == NULL ) buf = AFMALL(char, LBUF) ;
02190 
02191    /** step 1: read in ALL lines, see how many numbers are in each,
02192                in order to get the maximum row length and # of rows **/
02193 
02194    lbfill = filler ; /* 10 Aug 2004 */
02195 
02196    (void) my_fgets( NULL , 0 , NULL ) ;  /* reset */
02197    ncol = nrow = 0 ;
02198    while(1){
02199      ptr = my_fgets( buf , LBUF , fts ) ;
02200      if( ptr==NULL || *ptr=='\0' ) break ;
02201      fvec = decode_linebuf( buf ) ;
02202      if( fvec != NULL && fvec->nar > 0 ){ nrow++; ncol = MAX(ncol,fvec->nar); }
02203      if( fvec != NULL ) KILL_floatvec(fvec) ; else break ;
02204    }
02205    if( nrow == 0 || ncol == 0 ){ fclose(fts); FRB(buf); lbfill=0.0f; RETURN(NULL); }
02206 
02207    /** At this point, ncol is the number of floats to be read from each line **/
02208 
02209    rewind( fts ) ;  /* will start over */
02210 
02211    outim = mri_new( ncol , nrow , MRI_float ) ;
02212    tsar  = MRI_FLOAT_PTR(outim) ;
02213 
02214    /** read lines, convert to floats, store **/
02215 
02216    nrow = 0 ;
02217    while( 1 ){
02218      ptr = my_fgets( buf , LBUF , fts ) ;  /* read */
02219      if( ptr==NULL || *ptr=='\0' ) break ; /* failure --> end of data */
02220 
02221      fvec = decode_linebuf( buf ) ;
02222      if( fvec == NULL ) break ;
02223      if( fvec->nar == 0 ){ KILL_floatvec(fvec); break; }
02224 
02225      for( ii=0 ; ii < fvec->nar && ii < ncol ; ii++ )
02226        tsar[nrow*ncol+ii] = fvec->ar[ii] ;
02227      for( ; ii < ncol ; ii++ )
02228        tsar[nrow*ncol+ii] = filler ;   /* fill for incomplete lines */
02229      KILL_floatvec(fvec) ;
02230      nrow++ ;                  /* got one more complete row! */
02231    }
02232    fclose( fts ) ; /* finished with this file! */
02233    (void) my_fgets( NULL , 0 , NULL ) ;  /* reset */
02234 
02235    mri_add_name( fname , outim ) ;
02236    FRB(buf) ; lbfill = 0.0f ; RETURN(outim) ;
02237 }
02238 
02239 /*---------------------------------------------------------------------------
02240   Read in an ASCII file to a float array.
02241 -----------------------------------------------------------------------------*/
02242 
02243 static void read_ascii_floats( char * fname, int * nff , float ** ff )
02244 {
02245    int ii,jj,val , used_tsar , alloc_tsar ;
02246    float *tsar ;
02247    float ftemp ;
02248    FILE *fts ;
02249    char *buf ;  /* 08 Jul 2004: malloc this now, instead of auto */
02250    char *ptr ;
02251    int  bpos , blen , nrow ;
02252 
02253    /* check inputs */
02254 
02255    if( nff == NULL || ff == NULL ) return ;
02256    if( fname == NULL || fname[0] == '\0' ){ *nff=0 ; *ff=NULL ; return ; }
02257 
02258    fts = fopen( fname , "r" ) ;
02259    if( fts == NULL ){ *nff=0 ; *ff=NULL ; return ; }
02260 
02261    /* make some space */
02262 
02263    used_tsar  = 0 ;
02264    alloc_tsar = INC_TSARSIZE ;
02265    tsar       = (float *) malloc( sizeof(float) * alloc_tsar ) ;
02266    if( tsar == NULL ){
02267      fprintf(stderr,"\n*** malloc fails: read_ascii_floats ***\n"); EXIT(1);
02268    }
02269 
02270    /** read lines, convert to floats, store **/
02271 
02272    nrow = 0 ;
02273    buf = (char *)malloc(LBUF) ;
02274    while( 1 ){
02275       ptr = fgets( buf , LBUF , fts ) ;  /* read */
02276       if( ptr == NULL ) break ;          /* failure --> end of data */
02277       blen = strlen(buf) ;
02278       if( blen <= 0 ) break ;            /* nothing --> end of data */
02279 
02280       for( ii=0 ; ii < blen && isspace(buf[ii]) ; ii++ ) ; /* skip blanks */
02281 
02282       if( ii      == blen ) continue ;    /* skip all blank line */
02283       if( buf[ii] == '#'  ) continue ;    /* skip a comment line */
02284       if( buf[ii] == '!'  ) continue ;
02285 
02286       /* convert commas to blanks */
02287 
02288       for( jj=ii ; jj < blen ; jj++ ) if( buf[jj] == ',' ) buf[jj] = ' ' ;
02289 
02290       for( bpos=ii ; bpos < blen ; ){
02291          val = sscanf( buf+bpos , "%f%n" , &ftemp , &jj ) ;  /* read from string */
02292          if( val < 1 ) break ;                               /* bad read? */
02293          bpos += jj ;                                        /* start of next read */
02294 
02295          if( used_tsar == alloc_tsar ){
02296             alloc_tsar += INC_TSARSIZE ;
02297             tsar        = (float *)realloc( tsar,sizeof(float)*alloc_tsar );
02298             if( tsar == NULL ){
02299                fprintf(stderr,"\n*** realloc fails: read_ascii_floats ***\n"); EXIT(1);
02300             }
02301          }
02302 
02303          tsar[used_tsar++] = ftemp ;  /* store input */
02304       }
02305 
02306       nrow++ ;                  /* got one more complete row! */
02307    }
02308    fclose( fts ) ; /* finished with this file! */
02309    free( buf ) ;
02310 
02311    if( used_tsar <= 1 ){ free(tsar); *nff=0; *ff=NULL; return; }
02312 
02313    tsar = (float *) realloc( tsar , sizeof(float) * used_tsar ) ;
02314    if( tsar == NULL ){
02315       fprintf(stderr,"\n*** final realloc fails: read_ascii_floats ***\n"); EXIT(1);
02316    }
02317 
02318    *nff = used_tsar; *ff  = tsar; return;
02319 }
02320 
02321 /*--------------------------------------------------------------
02322    Read a pile of images from one ASCII file.
02323    Adapted from mri_read_3D - Feb 2000 - RWCox.
02324    [N.B.: if this routine is altered, don't forget mri_imcount!]
02325 ----------------------------------------------------------------*/
02326 
02327 MRI_IMARR * mri_read_3A( char * tname )
02328 {
02329    int nx , ny , nz , ii , nxyz,nxy , nff ;
02330    int ngood , length , kim , datum_type ;
02331    char fname[256]="\0" , buf[512] ;
02332    MRI_IMARR * newar ;
02333    MRI_IMAGE * newim , * flim ;
02334    float * ff ;
02335 
02336 ENTRY("mri_read_3A") ;
02337 
02338    /*** get info from 3A tname ***/
02339 
02340    if( tname == NULL || strlen(tname) < 10 ) RETURN(NULL) ;
02341 
02342    switch( tname[2] ){  /* allow for 3As:, 3Ab:, 3Af: */
02343 
02344       default: ngood = 0 ; break ;
02345 
02346       case 's':
02347          ngood = sscanf( tname, "3As:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
02348          datum_type = MRI_short ;
02349          break ;
02350 
02351       case 'b':
02352          ngood = sscanf( tname, "3Ab:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
02353          datum_type = MRI_byte ;
02354          break ;
02355 
02356       case 'f':
02357          ngood = sscanf( tname, "3Af:%d:%d:%d:%s", &nx, &ny, &nz, fname ) ;
02358          datum_type = MRI_float ;
02359          break ;
02360    }
02361 
02362    if( ngood < 4 || nx <= 0 || ny <= 0 || nz <= 0 || strlen(fname) <= 0 ) RETURN(NULL) ;
02363 
02364    /* read the input file */
02365 
02366    read_ascii_floats( fname , &nff , &ff ) ;
02367 
02368    if( nff <= 0 || ff == NULL ) RETURN(NULL) ;
02369 
02370    nxy = nx*ny ; nxyz = nxy*nz ;
02371 
02372    if( nff < nxyz ){
02373       fprintf(stderr,
02374                 "\n** WARNING: %s is too short - padding with %d zeros\n",
02375                 tname,nxyz-nff) ;
02376       ff = (float *) realloc( ff , sizeof(float) * nxyz ) ;
02377       for( ii=nff ; ii < nxyz ; ii++ ) ff[ii] = 0.0 ;
02378       nff = nxyz ;
02379    } else if( nff > nxyz ){
02380       fprintf(stderr,
02381                 "\n** WARNING: %s is too long - truncating off last %d values\n",
02382                 tname,nff-nxyz) ;
02383    }
02384 
02385    /* put the input data into MRI_IMAGEs */
02386 
02387    INIT_IMARR(newar) ;
02388 
02389    for( kim=0 ; kim < nz ; kim++ ){
02390       flim = mri_new( nx,ny , MRI_float ) ;
02391       memcpy( MRI_FLOAT_PTR(flim) , ff+nxy*kim , sizeof(float)*nxy ) ;
02392       switch( datum_type ){
02393          case MRI_float: newim = flim                                           ; break ;
02394          case MRI_short: newim = mri_to_short(1.0,flim)        ; mri_free(flim) ; break ;
02395          case MRI_byte:  newim = mri_to_byte_scl(1.0,0.0,flim) ; mri_free(flim) ; break ;
02396       }
02397 
02398       if( nz == 1 ) mri_add_name( fname , newim ) ;
02399       else {
02400          sprintf( buf , "%s#%d" , fname,kim ) ;
02401          mri_add_name( buf , newim ) ;
02402       }
02403 
02404       ADDTO_IMARR(newar,newim) ;
02405    }
02406 
02407    free(ff) ; RETURN(newar) ;
02408 }
02409 
02410 /*---------------------------------------------------------------------------
02411    Stuff to read an ANALYZE 7.5 .img file, given the .hdr filename
02412    -- 05 Feb 2001 - RWCox
02413    -- 27 Nov 2001 - modified to use funused1 as a scale factor
02414 -----------------------------------------------------------------------------*/
02415 
02416 #include "mayo_analyze.h"
02417 
02418 /*---------------------------------------------------------------*/
02419 /*! Byte swap ANALYZE file header in various places */
02420 
02421 static void swap_analyze_hdr( struct dsr *pntr )
02422 {
02423 ENTRY("swap_analyze_hdr") ;
02424    swap_4(&pntr->hk.sizeof_hdr) ;
02425    swap_4(&pntr->hk.extents) ;
02426    swap_2(&pntr->hk.session_error) ;
02427    swap_2(&pntr->dime.dim[0]) ;
02428    swap_2(&pntr->dime.dim[1]) ;
02429    swap_2(&pntr->dime.dim[2]) ;
02430    swap_2(&pntr->dime.dim[3]) ;
02431    swap_2(&pntr->dime.dim[4]) ;
02432    swap_2(&pntr->dime.dim[5]) ;
02433    swap_2(&pntr->dime.dim[6]) ;
02434    swap_2(&pntr->dime.dim[7]) ;
02435 #if 0
02436    swap_2(&pntr->dime.unused1) ;
02437 #endif
02438    swap_2(&pntr->dime.datatype) ;
02439    swap_2(&pntr->dime.bitpix) ;
02440    swap_4(&pntr->dime.pixdim[0]) ;
02441    swap_4(&pntr->dime.pixdim[1]) ;
02442    swap_4(&pntr->dime.pixdim[2]) ;
02443    swap_4(&pntr->dime.pixdim[3]) ;
02444    swap_4(&pntr->dime.pixdim[4]) ;
02445    swap_4(&pntr->dime.pixdim[5]) ;
02446    swap_4(&pntr->dime.pixdim[6]) ;
02447    swap_4(&pntr->dime.pixdim[7]) ;
02448    swap_4(&pntr->dime.vox_offset) ;
02449    swap_4(&pntr->dime.funused1) ;
02450    swap_4(&pntr->dime.funused2) ;
02451    swap_4(&pntr->dime.cal_max) ;
02452    swap_4(&pntr->dime.cal_min) ;
02453    swap_4(&pntr->dime.compressed) ;
02454    swap_4(&pntr->dime.verified) ;
02455    swap_2(&pntr->dime.dim_un0) ;
02456    swap_4(&pntr->dime.glmax) ;
02457    swap_4(&pntr->dime.glmin) ;
02458    EXRETURN ;
02459 }
02460 
02461 /*---------------------------------------------------------------*/
02462 /*! Count how many 2D slices are in an ANALYZE file.
02463 
02464    \param hname = the "hdr" file of the hdr/img file pair.
02465 */
02466 
02467 static int mri_imcount_analyze75( char * hname )
02468 {
02469    FILE * fp ;
02470    struct dsr hdr ;    /* ANALYZE .hdr format */
02471    int doswap , nz ;
02472 
02473 ENTRY("mri_imcount_analyze75") ;
02474 
02475    fp = fopen( hname , "rb" ) ;
02476    if( fp == NULL ) RETURN(0) ;
02477    hdr.dime.dim[0] = 0 ;
02478    fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
02479    fclose(fp) ;
02480    if( hdr.dime.dim[0] == 0 ) RETURN(0) ;
02481    doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
02482    if( doswap ) swap_analyze_hdr( &hdr ) ;
02483 
02484    switch( hdr.dime.dim[0] ){
02485       case 2:  nz = 1                                 ; break ;
02486       case 3:  nz = hdr.dime.dim[3]                   ; break ;
02487 
02488       default:
02489       case 4:  nz = hdr.dime.dim[3] * hdr.dime.dim[4] ; break ;
02490    }
02491    if( nz < 1 ) nz = 1 ;
02492 
02493    RETURN(nz) ;
02494 }
02495 
02496 /*---------------------------------------------------------------*/
02497 /*! Read an ANALYZE file into an ARRAY of 2D images.
02498 
02499    \param hname = the "hdr" file for the hdr/img pair
02500 */
02501 
02502 MRI_IMARR * mri_read_analyze75( char * hname )
02503 {
02504    FILE * fp ;
02505    char iname[1024] , buf[1024] ;
02506    int ii , jj , doswap ;
02507    struct dsr hdr ;    /* ANALYZE .hdr format */
02508    int ngood , length , kim , koff , datum_type , datum_len , swap ;
02509    int   nx,ny,nz , hglobal=0 , himage=0 ;
02510    float dx,dy,dz ;
02511    MRI_IMARR * newar ;
02512    MRI_IMAGE * newim ;
02513    void      * imar ;
02514    float fac=0.0 ;    /* 27 Nov 2001 */
02515    int floatize ;     /* 28 Nov 2001 */
02516    int spmorg=0 ;     /* 28 Nov 2001 */
02517 
02518 ENTRY("mri_read_analyze75") ;
02519 
02520    /* check & prepare filenames */
02521 
02522    if( hname == NULL ) RETURN(NULL) ;
02523    jj = strlen(hname) ;
02524    if( jj < 5 ) RETURN(NULL) ;
02525    if( strcmp(hname+jj-3,"hdr") != 0 ) RETURN(NULL) ;
02526    strcpy(iname,hname) ; strcpy(iname+jj-3,"img") ;
02527 
02528    /* read header file into struct */
02529 
02530    fp = fopen( hname , "rb" ) ;
02531    if( fp == NULL ) RETURN(NULL) ;
02532    hdr.dime.dim[0] = 0 ;
02533    fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
02534    fclose(fp) ;
02535    if( hdr.dime.dim[0] == 0 ) RETURN(NULL) ;
02536 
02537    /* check for swap-age */
02538 
02539    doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
02540    if( doswap ) swap_analyze_hdr( &hdr ) ;
02541 
02542    /* 28 Nov 2001: attempt to decode originator a la SPM */
02543 
02544    { short xyzuv[5] , xx,yy,zz ;
02545      memcpy( xyzuv , hdr.hist.originator , 10 ) ;
02546      if( xyzuv[3] == 0 && xyzuv[4] == 0 ){
02547         xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
02548         if( doswap ){ swap_2(&xx); swap_2(&yy); swap_2(&zz); }
02549         if( xx > 0 && xx < hdr.dime.dim[1] &&
02550             yy > 0 && yy < hdr.dime.dim[2] &&
02551             zz > 0 && zz < hdr.dime.dim[3]   ) spmorg = 1 ;
02552      }
02553    }
02554    if( spmorg ) strcpy( MRILIB_orients , "LRPAIS" ) ;
02555 
02556    /* 27 Nov 2001: get a scale factor for images */
02557 
02558    if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
02559       fac = hdr.dime.funused1 ;
02560       (void) thd_floatscan( 1 , &fac ) ;
02561       if( fac < 0.0 || fac == 1.0 ) fac = 0.0 ;
02562    }
02563 
02564    floatize = (fac != 0.0 || AFNI_yesenv("AFNI_ANALYZE_FLOATIZE")) ; /* 28 Nov 2001 */
02565 
02566    /* get data type into mrilib MRI_* form */
02567 
02568    switch( hdr.dime.datatype ){
02569       default:
02570          fprintf(stderr,"*** %s: Unknown ANALYZE datatype=%d (%s)\n",
02571                  hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
02572       RETURN(NULL) ;
02573 
02574       case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte   ;               break;
02575       case ANDT_SIGNED_SHORT:  datum_type = MRI_short  ;               break;
02576       case ANDT_SIGNED_INT:    datum_type = MRI_int    ;               break;
02577       case ANDT_FLOAT:         datum_type = MRI_float  ; floatize = 0; break;
02578       case ANDT_COMPLEX:       datum_type = MRI_complex; floatize = 0; break;
02579       case ANDT_RGB:           datum_type = MRI_rgb    ; floatize = 0; break;
02580       case ANDT_DOUBLE:        datum_type = MRI_double ; floatize = 0; break;
02581    }
02582 
02583    datum_len = mri_datum_size(datum_type) ;
02584 
02585    /* compute dimensions of images, and number of images */
02586 
02587    nx = hdr.dime.dim[1] ;
02588    ny = hdr.dime.dim[2] ;
02589    if( nx < 2 || ny < 2 ) RETURN(NULL) ;
02590 
02591    switch( hdr.dime.dim[0] ){
02592       case 2:  nz = 1                                 ; break ;
02593       case 3:  nz = hdr.dime.dim[3]                   ; break ;
02594 
02595       default:
02596       case 4:  nz = hdr.dime.dim[3] * hdr.dime.dim[4] ; break ;
02597    }
02598    if( nz < 1 ) nz = 1 ;
02599 
02600    dx = hdr.dime.pixdim[1] ;
02601    dy = hdr.dime.pixdim[2] ;
02602    dz = hdr.dime.pixdim[3] ;
02603 
02604    /** fprintf(stderr,"mri_read_analyze75: nx=%d ny=%d nz=%d\n",nx,ny,nz) ; **/
02605    /** fprintf(stderr,"mri_read_analyze75: dx=%g dy=%g dz=%g\n",dx,dy,dz) ; **/
02606 
02607    /* open .img file and read images from it */
02608 
02609    length = THD_filesize(iname) ;
02610    if( length <= 0 ){
02611       fprintf(stderr,"*** Can't find ANALYZE file %s\n",iname) ;
02612       RETURN(NULL) ;
02613    }
02614 
02615    fp = fopen( iname , "rb" ) ;
02616    if( fp == NULL ){
02617       fprintf(stderr,"*** Can't open ANALYZE file %s\n",iname) ;
02618       RETURN(NULL) ;
02619    }
02620 
02621    ngood = datum_len*nx*ny*nz ;
02622    if( length < ngood ){
02623       fprintf( stderr,
02624         "*** ANALYZE file %s is %d bytes long but must be at least %d bytes long\n"
02625         "*** for nx=%d ny=%d nz=%d and voxel=%d bytes\n",
02626         iname,length,ngood,nx,ny,nz,datum_len ) ;
02627       fclose(fp) ; RETURN(NULL) ;
02628    }
02629 
02630    /*** read images from the file ***/
02631 
02632    INIT_IMARR(newar) ;
02633 
02634    for( kim=0 ; kim < nz ; kim++ ){
02635       koff = hglobal + (kim+1)*himage + datum_len*nx*ny*kim ;
02636    /** fprintf(stderr,"mri_read_analyze75: kim=%d koff=%d\n",kim,koff) ; **/
02637       fseek( fp , koff , SEEK_SET ) ;
02638 
02639       newim  = mri_new( nx , ny , datum_type ) ;
02640       imar   = mri_data_pointer( newim ) ;
02641       length = fread( imar , datum_len , nx * ny , fp ) ;
02642 
02643       if( doswap ){
02644         switch( datum_len ){
02645           default: break ;
02646           case 2:  swap_twobytes (   nx*ny , imar ) ; break ;  /* short */
02647           case 4:  swap_fourbytes(   nx*ny , imar ) ; break ;  /* int, float */
02648           case 8:  swap_fourbytes( 2*nx*ny , imar ) ; break ;  /* complex */
02649         }
02650         newim->was_swapped = 1 ;  /* 07 Mar 2002 */
02651       }
02652 
02653       /* 28 Nov 2001: convert to floats? */
02654 
02655       if( floatize ){
02656          MRI_IMAGE *qim = mri_to_float(newim) ;
02657          mri_free(newim) ; newim = qim ;
02658       }
02659 
02660       if( nz == 1 ) mri_add_name( iname , newim ) ;
02661       else {
02662          sprintf( buf , "%s#%d" , iname,kim ) ;
02663          mri_add_name( buf , newim ) ;
02664       }
02665 
02666       newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dw = 1.0 ;
02667       ADDTO_IMARR(newar,newim) ;
02668 
02669       /* 27 Nov 2001: scale image? */
02670 
02671       if( fac != 0.0 ) mri_scale_inplace( fac , newim ) ;
02672    }
02673 
02674    fclose(fp) ; RETURN(newar) ;
02675 }
02676 
02677 /*-----------------------------------------------------------------*/
02678 /*! Read an ANALYZE file into an ARRAY of 3D images [26 Aug 2002].
02679 
02680    \param hname = the "hdr" file for the hdr/img pair
02681 */
02682 
02683 MRI_IMARR * mri_read3D_analyze75( char * hname )
02684 {
02685    FILE * fp ;
02686    char iname[1024] , buf[1024] ;
02687    int ii , jj , doswap ;
02688    struct dsr hdr ;    /* ANALYZE .hdr format */
02689    int ngood , length , kim , koff , datum_type , datum_len , swap ;
02690    int   nx,ny,nz , hglobal=0 , himage=0 ;
02691    float dx,dy,dz ;
02692    MRI_IMARR * newar ;
02693    MRI_IMAGE * newim ;
02694    void      * imar ;
02695    float fac=0.0 ;    /* 27 Nov 2001 */
02696    int floatize ;     /* 28 Nov 2001 */
02697    int spmorg=0 ;     /* 28 Nov 2001 */
02698 
02699    int   nt , nxyz ;  /* 26 Aug 2002 */
02700    float dt ;
02701 
02702 ENTRY("mri_read3D_analyze75") ;
02703 
02704    /* check & prepare filenames */
02705 
02706    if( hname == NULL ) RETURN(NULL) ;
02707    jj = strlen(hname) ;
02708    if( jj < 5 ) RETURN(NULL) ;
02709    if( strcmp(hname+jj-3,"hdr") != 0 ) RETURN(NULL) ;
02710    strcpy(iname,hname) ; strcpy(iname+jj-3,"img") ;
02711 
02712    /* read header file into struct */
02713 
02714    fp = fopen( hname , "rb" ) ;
02715    if( fp == NULL ) RETURN(NULL) ;
02716    hdr.dime.dim[0] = 0 ;
02717    fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
02718    fclose(fp) ;
02719    if( hdr.dime.dim[0] == 0 ) RETURN(NULL) ;
02720 
02721    /* check for swap-age */
02722 
02723    doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
02724    if( doswap ) swap_analyze_hdr( &hdr ) ;
02725 
02726    /* 28 Nov 2001: attempt to decode originator a la SPM */
02727 
02728    { short xyzuv[5] , xx,yy,zz ;
02729      memcpy( xyzuv , hdr.hist.originator , 10 ) ;
02730      if( xyzuv[3] == 0 && xyzuv[4] == 0 ){
02731         xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
02732         if( doswap ){ swap_2(&xx); swap_2(&yy); swap_2(&zz); }
02733         if( xx > 0 && xx < hdr.dime.dim[1] &&
02734             yy > 0 && yy < hdr.dime.dim[2] &&
02735             zz > 0 && zz < hdr.dime.dim[3]   ) spmorg = 1 ;
02736      }
02737    }
02738    if( spmorg ) strcpy( MRILIB_orients , "LRPAIS" ) ;
02739 
02740    /* 27 Nov 2001: get a scale factor for images */
02741 
02742    if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
02743       fac = hdr.dime.funused1 ;
02744       (void) thd_floatscan( 1 , &fac ) ;
02745       if( fac < 0.0 || fac == 1.0 ) fac = 0.0 ;
02746    }
02747 
02748    floatize = (fac != 0.0 || AFNI_yesenv("AFNI_ANALYZE_FLOATIZE")) ; /* 28 Nov 2001 */
02749 
02750    /* get data type into mrilib MRI_* form */
02751 
02752    switch( hdr.dime.datatype ){
02753       default:
02754          fprintf(stderr,"*** %s: Unknown ANALYZE datatype=%d (%s)\n",
02755                  hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
02756       RETURN(NULL) ;
02757 
02758       case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte   ;               break;
02759       case ANDT_SIGNED_SHORT:  datum_type = MRI_short  ;               break;
02760       case ANDT_SIGNED_INT:    datum_type = MRI_int    ;               break;
02761       case ANDT_FLOAT:         datum_type = MRI_float  ; floatize = 0; break;
02762       case ANDT_COMPLEX:       datum_type = MRI_complex; floatize = 0; break;
02763       case ANDT_RGB:           datum_type = MRI_rgb    ; floatize = 0; break;
02764    }
02765 
02766    datum_len = mri_datum_size(datum_type) ;
02767 
02768    /* compute dimensions of images, and number of images */
02769 
02770    nx = hdr.dime.dim[1] ;
02771    ny = hdr.dime.dim[2] ;
02772    if( nx < 2 || ny < 2 ) RETURN(NULL) ;
02773 
02774    switch( hdr.dime.dim[0] ){
02775       case 2:  nz = 1 ; nt = 1 ;                           ; break ;
02776       case 3:  nz = hdr.dime.dim[3] ; nt = 1 ;             ; break ;
02777 
02778       default:
02779       case 4:  nz = hdr.dime.dim[3] ; nt = hdr.dime.dim[4] ; break ;
02780    }
02781    if( nz < 1 ) nz = 1 ;
02782    if( nt < 1 ) nt = 1 ;
02783 
02784    dx = hdr.dime.pixdim[1] ;
02785    dy = hdr.dime.pixdim[2] ;
02786    dz = hdr.dime.pixdim[3] ;
02787    dt = hdr.dime.pixdim[4] ; if( dt <= 0.0 ) dt = 1.0 ;
02788 
02789    /* open .img file and read images from it */
02790 
02791    length = THD_filesize(iname) ;
02792    if( length <= 0 ){
02793       fprintf(stderr,"*** Can't find ANALYZE file %s\n",iname) ;
02794       RETURN(NULL) ;
02795    }
02796 
02797    fp = fopen( iname , "rb" ) ;
02798    if( fp == NULL ){
02799       fprintf(stderr,"*** Can't open ANALYZE file %s\n",iname) ;
02800       RETURN(NULL) ;
02801    }
02802 
02803    ngood = datum_len*nx*ny*nz*nt ;
02804    if( length < ngood ){
02805       fprintf( stderr,
02806         "*** ANALYZE file %s is %d bytes long but must be at least %d bytes long\n"
02807         "*** for nx=%d ny=%d nz=%d nt=%d and voxel=%d bytes\n",
02808         iname,length,ngood,nx,ny,nz,nt,datum_len ) ;
02809       fclose(fp) ; RETURN(NULL) ;
02810    }
02811 
02812    /*** read images from the file ***/
02813 
02814    INIT_IMARR(newar) ;
02815 
02816    for( kim=0 ; kim < nt ; kim++ ){
02817       koff = hglobal + (kim+1)*himage + datum_len*nxyz*kim ;
02818       fseek( fp , koff , SEEK_SET ) ;
02819 
02820       newim  = mri_new_vol( nx,ny,nz , datum_type ) ;
02821       imar   = mri_data_pointer( newim ) ;
02822       length = fread( imar , datum_len , nxyz , fp ) ;
02823 
02824       if( doswap ){
02825         switch( datum_len ){
02826           default: break ;
02827           case 2:  swap_twobytes (   nxyz , imar ) ; break ;  /* short */
02828           case 4:  swap_fourbytes(   nxyz , imar ) ; break ;  /* int, float */
02829           case 8:  swap_fourbytes( 2*nxyz , imar ) ; break ;  /* complex */
02830         }
02831         newim->was_swapped = 1 ;  /* 07 Mar 2002 */
02832       }
02833 
02834       /* 28 Nov 2001: convert to floats? */
02835 
02836       if( floatize ){
02837          MRI_IMAGE *qim = mri_to_float(newim) ;
02838          mri_free(newim) ; newim = qim ;
02839       }
02840 
02841       if( nt == 1 ) mri_add_name( iname , newim ) ;
02842       else {
02843          sprintf( buf , "%s#%d" , iname,kim ) ;
02844          mri_add_name( buf , newim ) ;
02845       }
02846 
02847       newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dt = dt ; newim->dw = 1.0 ;
02848       ADDTO_IMARR(newar,newim) ;
02849 
02850       /* 27 Nov 2001: scale image? */
02851 
02852       if( fac != 0.0 ) mri_scale_inplace( fac , newim ) ;
02853    }
02854 
02855    fclose(fp) ; RETURN(newar) ;
02856 }
02857 
02858 /*---------------------------------------------------------------------------
02859   12 Mar 2001 - stuff to read a Siemens Vision .ima file
02860 -----------------------------------------------------------------------------*/
02861 
02862 #include "siemens_vision.h"
02863 
02864 /*! Count the number of 2D images in a Siemens Vision .ima file.
02865 
02866    Unfortunately, this requires reading the image data and checking
02867    for all-zero images.  This is because Siemens stores their data
02868    in a fixed size file, and so just fills out the empty space with
02869    blank images if need be.
02870 */
02871 
02872 static int mri_imcount_siemens( char * hname )
02873 {
02874    struct Siemens_vision_header head ;
02875    FILE * fp ;
02876    int i,j,xx,yy , matrix , swap , imagesize,nxx,blank , slices ;
02877    struct stat file_stat ;
02878    short *imar ;
02879 
02880    /*--- check file size ---*/
02881 
02882    if( hname == NULL ) return 0 ;
02883 
02884    i = stat( hname , &file_stat ) ;
02885    if( i < 0 ) return 0 ;
02886 
02887    /*--- read header data ---*/
02888 
02889    fp = fopen( hname , "rb" ) ;
02890    if( fp == NULL ) return 0 ;
02891    fread( &head , sizeof(struct Siemens_vision_header) , 1 , fp ) ;
02892 
02893    /*-- check some integer in header to determine if we need to byteswap --*/
02894 
02895    swap = ( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ) ;
02896    if( swap ){
02897       swap_4( &(head.SiemensStudyDateMM) ) ;
02898       if( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ){
02899          swap = 0 ;
02900       }
02901    }
02902 
02903    /*-- find image size from header --*/
02904 
02905    if( swap ) swap_4( &(head.DisplayMatrixSize) ) ;
02906    imagesize = head.DisplayMatrixSize ;
02907 
02908    /*-- determine number of sub-images in file --*/
02909 
02910 #undef  MATRIX_MAX
02911 #define MATRIX_MAX 9
02912 
02913    i = 2*imagesize*imagesize ;
02914    for( matrix=1 ; matrix < MATRIX_MAX ; matrix++ )
02915      if( file_stat.st_size == i*matrix*matrix + SIEMENS_HEADERSIZE ) break ;
02916 
02917    if( matrix == MATRIX_MAX ){
02918      fclose(fp) ; return 0 ; /* didn't recognize file format */
02919    }
02920 #undef MATRIX_MAX
02921 
02922    /*-- read image data from file (but don't byteswap it) --*/
02923 
02924    imar = (short *) calloc(sizeof(short),matrix*matrix*imagesize*imagesize) ;
02925    fseek( fp , SIEMENS_HEADERSIZE , SEEK_SET ) ;
02926    fread( imar , sizeof(short) , matrix*matrix*imagesize*imagesize , fp ) ;
02927    fclose(fp) ;
02928 
02929    /*-- count slices - all zero (blank) slices at end are skipped --*/
02930 
02931    slices = 0 ; nxx = matrix*imagesize ;
02932 
02933    for( yy=0 ; yy < matrix ; yy++ ){      /* rows in array of sub-images */
02934       for( xx=0 ; xx < matrix ; xx++ ){   /* cols in array of sub-images */
02935          blank = 1 ;
02936          for( j=0 ; j < imagesize ; j++ ){    /* row in sub-image */
02937             for( i=0 ; i < imagesize ; i++ ){ /* col in sub-image */
02938                if( imar[i+xx*imagesize+(j+yy*imagesize)*nxx] ) blank = 0 ;
02939             }
02940          }
02941          if( !blank ) slices = 1 + xx + yy*matrix ;
02942       }
02943    }
02944 
02945    free(imar) ; return slices ;
02946 }
02947 
02948 /*---------------------------------------------------------------------------*/
02949 /*! Read an array of 2D images from Siemens Vision .ima file.
02950 
02951    The images are stored in a 2D array, which requires untangling the
02952    data rows to put them into separate MRI_IMAGE structs.
02953 */
02954 
02955 MRI_IMARR * mri_read_siemens( char * hname )
02956 {
02957    struct Siemens_vision_header head ;
02958    FILE * fp ;
02959    int i,j,xx,yy , matrix , swap , imagesize,nxx,blank , slices,nz ;
02960    struct stat file_stat ;
02961    short *imar ;
02962    MRI_IMARR * newar ;
02963    MRI_IMAGE * newim ;
02964    short     * nar ;
02965    char buf[256] ;
02966    float dx,dy,dz ;
02967    char *eee ; int ileave=0 ;  /* 25 Sep 2001 */
02968 
02969 ENTRY("mri_read_siemens") ;
02970 
02971    /*--- check file size ---*/
02972 
02973    if( hname == NULL ) RETURN(NULL) ;
02974 
02975    i = stat( hname , &file_stat ) ;
02976    if( i < 0 ) RETURN(NULL) ;
02977 
02978    /*--- read header data ---*/
02979 
02980    fp = fopen( hname , "rb" ) ;
02981    if( fp == NULL ) RETURN(NULL) ;
02982    fread( &head , sizeof(struct Siemens_vision_header) , 1 , fp ) ;
02983 
02984    /*-- check some integer in header to determine if we need to byteswap --*/
02985 
02986    swap = ( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ) ;
02987    if( swap ){
02988       swap_4( &(head.SiemensStudyDateMM) ) ;
02989       if( head.SiemensStudyDateMM < 0 || head.SiemensStudyDateMM > 13 ){
02990          swap = 0 ;
02991       }
02992    }
02993 
02994    /*-- find image size from header --*/
02995 
02996    if( swap ) swap_4( &(head.DisplayMatrixSize) ) ;
02997    imagesize = head.DisplayMatrixSize ;
02998 
02999    /*-- determine number of sub-images in file --*/
03000 
03001 #undef  MATRIX_MAX
03002 #define MATRIX_MAX 16
03003 
03004    i = 2*imagesize*imagesize ;
03005    for( matrix=1 ; matrix < MATRIX_MAX ; matrix++ )
03006      if( file_stat.st_size == i*matrix*matrix + SIEMENS_HEADERSIZE ) break ;
03007 
03008    if( matrix == MATRIX_MAX ){
03009      fclose(fp) ; RETURN(NULL) ; /* didn't recognize file format */
03010    }
03011 #undef MATRIX_MAX
03012 
03013    /*-- read image data from file and byteswap it, if needed --*/
03014 
03015    imar = (short *) calloc(sizeof(short),matrix*matrix*imagesize*imagesize) ;
03016    fseek( fp , SIEMENS_HEADERSIZE , SEEK_SET ) ;
03017    fread( imar , sizeof(short) , matrix*matrix*imagesize*imagesize , fp ) ;
03018    fclose(fp) ;
03019 
03020    if( swap ) swap_twobytes( matrix*matrix*imagesize*imagesize , imar ) ;
03021 
03022    /*-- count slices - all zero (blank) slices at end are skipped --*/
03023 
03024    slices = 0 ; nxx = matrix*imagesize ;
03025 
03026    for( yy=0 ; yy < matrix ; yy++ ){      /* rows in array of sub-images */
03027       for( xx=0 ; xx < matrix ; xx++ ){   /* cols in array of sub-images */
03028          blank = 1 ;
03029          for( j=0 ; j < imagesize ; j++ ){    /* row in sub-image */
03030             for( i=0 ; i < imagesize ; i++ ){ /* col in sub-image */
03031                if( imar[i+xx*imagesize+(j+yy*imagesize)*nxx] ) blank = 0 ;
03032             }
03033          }
03034          if( !blank ) slices = 1 + xx + yy*matrix ;
03035       }
03036    }
03037 
03038    if( slices == 0 ){ free(imar) ; RETURN(NULL) ; }  /* bad news */
03039 
03040    /*-- get image dimensions, etc --*/
03041 
03042    if( swap ){
03043      swap_8(&(head.FOVRow));
03044      swap_8(&(head.FOVColumn));
03045      swap_8(&(head.SliceThickness));
03046    }
03047    dx = head.FOVRow    / imagesize ;
03048    dy = head.FOVColumn / imagesize ;
03049    dz = head.SliceThickness ;
03050 
03051    /*-- save orientation and offset in global variables --*/
03052 
03053    MRILIB_orients[0] = head.OrientationSet1Left[0] ;
03054    MRILIB_orients[1] = head.OrientationSet2Right[0];
03055    MRILIB_orients[2] = head.OrientationSet1Top[0]  ;
03056    MRILIB_orients[3] = head.OrientationSet2Down[0] ;
03057    MRILIB_orients[4] = head.OrientationSet1Back[0] ;
03058    MRILIB_orients[5] = head.OrientationSet2Front[0];
03059    for (i=0; i<6; i++) {
03060      if (MRILIB_orients[i]=='H') MRILIB_orients[i]='S';
03061      if (MRILIB_orients[i]=='F') MRILIB_orients[i]='I';
03062    }
03063    MRILIB_orients[6] = '\0' ;
03064    MRILIB_zoff = fabs(strtod(head.TextSlicePosition,NULL)) ; use_MRILIB_zoff = 1 ;
03065 
03066    /*-- create output --*/
03067 
03068    INIT_IMARR(newar) ;
03069 
03070    for( yy=0 ; yy < matrix ; yy++ ){      /* rows in array of sub-images */
03071       for( xx=0 ; xx < matrix ; xx++ ){   /* cols in array of sub-images */
03072 
03073          newim = mri_new( imagesize , imagesize , MRI_short ) ;
03074          nar   = MRI_SHORT_PTR( newim ) ;
03075 
03076          if( swap ) newim->was_swapped = 1 ; /* 07 Mar 2002 */
03077 
03078          for( j=0 ; j < imagesize ; j++ )    /* row in sub-image */
03079            memcpy( nar+j*imagesize ,
03080                    imar+xx*imagesize+(j+yy*imagesize)*nxx , 2*imagesize ) ;
03081 
03082          sprintf( buf , "%s#%d:%d" , hname,xx,yy ) ;
03083          mri_add_name( buf , newim ) ;
03084 
03085          newim->dx = dx ; newim->dy = dy ; newim->dz = dz ; newim->dw = 1.0 ;
03086          ADDTO_IMARR(newar,newim) ;
03087          if( IMARR_COUNT(newar) == slices ) goto Done ;  /* Aauugghh!!! */
03088       }
03089    }
03090 
03091 Done:
03092 
03093    /*-- 25 Sep 2001: possibly interleave the images --*/
03094 
03095    eee = getenv("AFNI_SIEMENS_INTERLEAVE") ;
03096    ileave = ( (eee != NULL) && (*eee=='Y' || *eee=='y') ) ;
03097    if( ileave && slices > 2 ){
03098       int mid = (slices-1)/2 ;  /* midpoint */
03099       MRI_IMARR *qar ;          /* new image array */
03100       INIT_IMARR(qar) ;
03101       for( i=0 ; i < slices ; i++ ){
03102          if( i%2 == 0 ) j = i/2 ;           /* slice #i is in newar #j */
03103          else           j = mid + (i+1)/2 ;
03104          ADDTO_IMARR(qar,IMARR_SUBIM(newar,j)) ; /* move image to new array */
03105       }
03106       FREE_IMARR(newar) ; newar = qar ;
03107    }
03108 
03109    free(imar) ; RETURN(newar) ;
03110 }
03111 
03112 /*---------------------------------------------------------------------------*/
03113 /*! Check for dicom magic number (string) in file
03114     Bytes 128-131 should be "DICM" in a Dicom Part 10 file
03115 */
03116 
03117 short check_dicom_magic_num( char * fname )
03118 {
03119   FILE * fp;
03120   char test_string[5] ;
03121 
03122   fp = fopen( fname, "rb" ) ;
03123   if(fp == NULL ) return 0 ;
03124   fseek( fp, 128 , SEEK_SET ) ;
03125   fread( test_string , 1 , 4 , fp ) ; test_string[4] = '\0' ;
03126   fclose( fp ) ;
03127   if( strcmp(test_string,"DICM") == 0 ) {
03128     return 1 ;
03129   } else {
03130     return 0 ;
03131   }
03132 }
03133 
03134 /*---------------------------------------------------------------------------
03135    Stuff to read a file in "delay" mode -- 01 Jan 1997.
03136 -----------------------------------------------------------------------------*/
03137 
03138 #ifdef USE_MRI_DELAY
03139 
03140 /**** If possible, throw the data away for later retrieval from disk ****/
03141 
03142 void mri_purge_delay( MRI_IMAGE * im )
03143 {
03144    void * ar ;
03145 
03146    /** if no delay filename,
03147        or if it is marked as already set for delay input, do nothing **/
03148 
03149    if( im->fname == NULL ||
03150        (im->fondisk & INPUT_DELAY) != 0 ) return ;
03151 
03152    /** get the data pointer, throw data way, clear the data pointer **/
03153 
03154    ar = mri_data_pointer( im ) ;
03155    if( ar != NULL ){ free(ar) ; mri_clear_data_pointer(im) ; }
03156 
03157    /** mark as set for delay input **/
03158 
03159    im->fondisk |= INPUT_DELAY ;
03160    return ;
03161 }
03162 
03163 /**** if possible, read data from delay input file ****/
03164 
03165 void mri_input_delay( MRI_IMAGE * im )
03166 {
03167    FILE * imfile=NULL ;
03168    void * imar ;
03169 
03170    /** if no delay input file,
03171        or is marked as already read in, do nothing **/
03172 
03173    if( im->fname == NULL ||
03174        (im->fondisk & INPUT_DELAY) == 0 ) return ;
03175 
03176    /** open the delay input file [06 Mar 2001: maybe not] **/
03177 
03178    if( strcmp(im->fname,"ALLZERO") != 0 ){
03179       imfile = fopen( im->fname , "r" ) ;
03180       if( imfile == NULL ){
03181          fprintf( stderr , "couldn't open delayed image file %s\n" , im->fname ) ;
03182          return ;
03183       }
03184    }
03185 
03186    /** make space for the array **/
03187 
03188    imar = (void *) malloc( im->nvox * im->pixel_size ) ;
03189    if( imar == NULL ){
03190       fprintf( stderr ,
03191                "malloc fails for delayed image from file %s\n" , im->fname ) ;
03192       if( imfile != NULL ) fclose( imfile ) ;
03193       return ;
03194    }
03195    mri_fix_data_pointer( imar , im ) ;
03196 
03197    /** read from the file into the array **/
03198 
03199    if( imfile != NULL ){
03200       fseek( imfile , im->foffset , SEEK_SET ) ;
03201       fread( imar , im->pixel_size , im->nvox , imfile ) ;
03202       fclose( imfile ) ;
03203    } else {
03204       memset( imar , 0 , im->nvox * im->pixel_size ) ;  /* 06 Mar 2001 */
03205    }
03206 
03207    /** swap bytes, if so marked **/
03208 
03209    if( (im->fondisk & BSWAP_DELAY) ){
03210       mri_swapbytes( im ) ;
03211       im->was_swapped = 1 ;  /* 07 Mar 2002 */
03212    }
03213 
03214    /** mark as already read from disk **/
03215 
03216    im->fondisk ^= INPUT_DELAY ;
03217 
03218 #if 0
03219 fprintf(stderr,"delayed input from file %s at offset %d\n",im->fname,im->foffset);
03220 #endif
03221    return ;
03222 }
03223 
03224 /**********************************************************************/
03225 /**** like mri_read_file, but returns delayed images for 3D: files ****/
03226 /**** (all others are read in now anyhoo, so there)                ****/
03227 
03228 MRI_IMARR * mri_read_file_delay( char * fname )
03229 {
03230    MRI_IMARR * newar=NULL ;
03231    MRI_IMAGE * newim ;
03232    char * new_fname ;
03233 
03234    new_fname = imsized_fname( fname ) ;
03235    if( new_fname == NULL ) return NULL ;
03236 
03237    if( strlen(new_fname) > 9 && new_fname[0] == '3' && new_fname[1] == 'D' &&
03238        (new_fname[2] == ':' || new_fname[3] == ':') ){
03239                                /* check for ':', too   3 Jan 2005 [rickr] */
03240 
03241       newar = mri_read_3D_delay( new_fname ) ;   /* read from a 3D file, later */
03242 
03243    } else if( strlen(new_fname) > 9 &&
03244               new_fname[0] == '3' && new_fname[1] == 'A' && new_fname[3] == ':' ){
03245 
03246       newar = mri_read_3A( new_fname ) ;
03247 
03248    } else if( strstr(new_fname,".hdr") != NULL ||
03249               strstr(new_fname,".HDR") != NULL   ){ /* 05 Feb 2001 - ANALYZE header */
03250 
03251       newar = mri_read_analyze75( new_fname ) ;
03252 
03253    } else if( strstr(new_fname,".ima") != NULL ||
03254               strstr(new_fname,".IMA") != NULL   ){ /* 12 Mar 2001 - Siemens */
03255 
03256       newar = mri_read_siemens( new_fname ) ;
03257 
03258    } else if( strstr(new_fname,".mpg" ) != NULL ||  /* 03 Dec 2003 */
03259               strstr(new_fname,".MPG" ) != NULL ||  /* read MPEGs  */
03260               strstr(new_fname,".mpeg") != NULL ||
03261               strstr(new_fname,".MPEG") != NULL   ){
03262 
03263       newar = mri_read_mpeg( new_fname ) ;  /* cf. mri_read_mpeg.c */
03264    }
03265 
03266    /* failed thus far?  try DICOM, unless user has requested DICOM last */
03267    /* 05 May 2003 added option to try DICOM last         KRH          */
03268 
03269    if ((newar == NULL) && !AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
03270      newar = mri_read_dicom( new_fname ) ;  /* cf. mri_read_dicom.c */
03271    }
03272 
03273    /* failed again?  try mri_read() for 1 image */
03274 
03275    if( newar == NULL ){
03276       newim = mri_read( new_fname ) ;      /* read from a 2D file */
03277       if( newim == NULL ){ free(new_fname) ; return NULL ; }
03278       INIT_IMARR(newar) ;
03279       ADDTO_IMARR(newar,newim) ;
03280    }
03281 
03282    if ( (newar == NULL) && AFNI_yesenv("AFNI_TRY_DICOM_LAST")) {
03283      newar = mri_read_dicom( new_fname ) ;  /* cf. mri_read_dicom.c */
03284    }
03285 
03286    free(new_fname) ;
03287    return newar ;
03288 }
03289 
03290 /**** like mri_read_3D, but returns delayed images ****/
03291 
03292 MRI_IMARR * mri_read_3D_delay( char * tname )
03293 {
03294    int hglobal , himage , nx , ny , nz ;
03295    char fname[256] , buf[512] ;
03296    int ngood , length , kim , datum_type , datum_len , swap ;
03297    MRI_IMARR * newar ;
03298    MRI_IMAGE * newim ;
03299    FILE      * imfile ;
03300 
03301    /*** get info from 3D tname ***/
03302 
03303    if( tname == NULL || strlen(tname) < 10 ) return NULL ;
03304 
03305    switch( tname[2] ){  /* allow for 3D: or 3Ds: or 3Db: */
03306 
03307       default:
03308       case ':':
03309          ngood = sscanf( tname , "3D:%d:%d:%d:%d:%d:%s" ,
03310                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03311 
03312          swap       = 0 ;
03313          datum_type = MRI_short ;
03314          datum_len  = sizeof(short) ;  /* better be 2 */
03315          break ;
03316 
03317       case 's':
03318          ngood = sscanf( tname , "3Ds:%d:%d:%d:%d:%d:%s" ,
03319                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03320 
03321          swap       = 1 ;
03322          datum_type = MRI_short ;
03323          datum_len  = sizeof(short) ;  /* better be 2 */
03324          break ;
03325 
03326       case 'b':
03327          ngood = sscanf( tname , "3Db:%d:%d:%d:%d:%d:%s" ,
03328                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03329 
03330          swap       = 0 ;
03331          datum_type = MRI_byte ;
03332          datum_len  = sizeof(byte) ;  /* better be 1 */
03333          break ;
03334 
03335       case 'f':
03336          ngood = sscanf( tname , "3Df:%d:%d:%d:%d:%d:%s" ,
03337                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03338 
03339          swap       = 0 ;
03340          datum_type = MRI_float ;
03341          datum_len  = sizeof(float) ;  /* better be 4 */
03342          break ;
03343 
03344       case 'd':                                            /* 06 Feb 2003 */
03345          ngood = sscanf( tname , "3Dd:%d:%d:%d:%d:%d:%s" ,
03346                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03347 
03348          swap       = 0 ;
03349          datum_type = MRI_float ;
03350          datum_len  = sizeof(double) ;  /* better be 8 */
03351          break ;
03352 
03353       case 'i':
03354          ngood = sscanf( tname , "3Di:%d:%d:%d:%d:%d:%s" ,
03355                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03356 
03357          swap       = 0 ;
03358          datum_type = MRI_int ;
03359          datum_len  = sizeof(int) ;  /* better be 4 */
03360          break ;
03361 
03362       case 'c':
03363          ngood = sscanf( tname , "3Dc:%d:%d:%d:%d:%d:%s" ,
03364                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03365 
03366          swap       = 0 ;
03367          datum_type = MRI_complex ;
03368          datum_len  = sizeof(complex) ;  /* better be 8 */
03369          break ;
03370 
03371       case 'r':
03372          ngood = sscanf( tname , "3Dr:%d:%d:%d:%d:%d:%s" ,
03373                          &hglobal , &himage , &nx , &ny , &nz , fname ) ;
03374 
03375          swap       = 0 ;
03376          datum_type = MRI_rgb ;
03377          datum_len  = 3*sizeof(byte) ;  /* better be 3 */
03378          break ;
03379    }
03380 
03381    if( ngood < 6 || himage < 0 ||
03382        nx <= 0   || ny <= 0    || nz <= 0 ||
03383        strlen(fname) <= 0                   ) return NULL ;   /* bad info */
03384 
03385    /*** open the input file and position it [06 Mar 2001: maybe not] ***/
03386 
03387    if( strcmp(fname,"ALLZERO") != 0 ){
03388       imfile = fopen( fname , "r" ) ;
03389       if( imfile == NULL ){
03390          fprintf( stderr , "couldn't open delayed image file %s\n" , fname ) ;
03391          return NULL ;
03392       }
03393    } else {
03394       imfile = NULL ;
03395    }
03396 
03397    if( imfile != NULL ){
03398       fseek( imfile , 0L , SEEK_END ) ;  /* get the length of the file */
03399       length = ftell( imfile ) ;
03400 
03401    /** 13 Apr 1999: modified to allow actual hglobal < -1
03402                     as long as hglobal+himage >= 0       **/
03403 
03404 #if 0                 /* the old code */
03405       if( hglobal < 0 ){
03406          hglobal = length - nz*(datum_len*nx*ny+himage) ;
03407          if( hglobal < 0 ) hglobal = 0 ;
03408       }
03409 #else                 /* 13 Apr 1999 */
03410       if( hglobal == -1 || hglobal+himage < 0 ){
03411          hglobal = length - nz*(datum_len*nx*ny+himage) ;
03412          if( hglobal < 0 ) hglobal = 0 ;
03413       }
03414 #endif
03415 
03416       ngood = hglobal + nz*(datum_len*nx*ny+himage) ;
03417       if( length < ngood ){
03418          fprintf( stderr,
03419            "file %s is %d bytes long but must be at least %d bytes long\n"
03420            "for hglobal=%d himage=%d nx=%d ny=%d nz=%d and voxel=%d bytes\n",
03421            fname,length,ngood,hglobal,himage,nx,ny,nz,datum_len ) ;
03422          fclose( imfile ) ;
03423          return NULL ;
03424       }
03425       fclose( imfile ) ;
03426    }
03427 
03428    /*** put pointers to data in the file into the images ***/
03429 
03430    INIT_IMARR(newar) ;
03431 
03432    for( kim=0 ; kim < nz ; kim++ ){
03433       newim = mri_new_vol_empty( nx,ny,1 , datum_type ) ;  /* empty image */
03434       mri_add_fname_delay( fname , newim ) ;               /* put filename in */
03435       newim->fondisk = (swap) ? (INPUT_DELAY | BSWAP_DELAY) /* mark read type */
03436                               : (INPUT_DELAY) ;
03437       newim->foffset = hglobal + (kim+1)*himage + datum_len*nx*ny*kim ;
03438 
03439       if( nz == 1 ) mri_add_name( fname , newim ) ;
03440       else {
03441          sprintf( buf , "%s#%d" , fname,kim ) ;
03442          mri_add_name( buf , newim ) ;
03443       }
03444 
03445       ADDTO_IMARR(newar,newim) ;
03446    }
03447 
03448    return newar ;
03449 }
03450 #endif /* USE_MRI_DELAY */