00001 #include "mrilib.h"
00002 
00003 /*-----------------------------------------------------------------------------------*/
00004 
00005 #define NMOMAX 256
00006 typedef struct {
00007   int good ;                       /* data in here is good? */
00008   int have_data[3] ;               /* do we have slices 0 and 1 in           *
00009                                     * each dimension to determine z-spacing? *
00010                                     * added 25 Feb 2003 KRH                  */
00011   int mosaic_num ;                 /* how many slices in 1 'image' */
00012   float slice_xyz[NMOMAX][3] ;     /* Sag, Cor, Tra coordinates */
00013 } Siemens_extra_info ;
00014 
00015 static char * extract_bytes_from_file( FILE *fp, off_t start, size_t len, int strize ) ;
00016 static void get_siemens_extra_info( char *str , Siemens_extra_info *mi ) ;
00017 
00018 /*-----------------------------------------------------------------------------------*/
00019 /* Save the Siemens extra info string in case the caller wants to get it. */
00020 
00021 static char *str_sexinfo=NULL ;
00022 
00023 char * mri_dicom_sexinfo(void){ return str_sexinfo; }  /* 23 Dec 2002 */
00024 
00025 /*-----------------------------------------------------------------------------------*/
00026 
00027 static int LITTLE_ENDIAN_ARCHITECTURE = -1 ;
00028 
00029 static void RWC_set_endianosity(void)
00030 {
00031    union { unsigned char bb[2] ;
00032            short         ss    ; } fred ;
00033 
00034    if( LITTLE_ENDIAN_ARCHITECTURE < 0 ){
00035      fred.bb[0] = 1 ; fred.bb[1] = 0 ;
00036 
00037      LITTLE_ENDIAN_ARCHITECTURE = (fred.ss == 1) ;
00038    }
00039 }
00040 
00041 /*-----------------------------------------------------------------------------------*/
00042 
00043 static char *elist[] = {
00044 
00045  "0018 0050" ,  /* Slice thickness */
00046  "0018 0080" ,  /* Repetition time */
00047  "0018 0088" ,  /* Spacing between slices */
00048  "0018 1149" ,  /* Field of view */
00049 
00050  "0020 0020" ,  /* Patient orientation */
00051  "0020 0032" ,  /* Image position (patient) */
00052  "0020 0037" ,  /* Image orientation (patient) */
00053  "0020 1041" ,  /* Slice location */
00054 
00055  "0028 0002" ,  /* Samples per pixel */
00056  "0028 0008" ,  /* Number of frames */
00057  "0028 0010" ,  /* Rows */
00058  "0028 0011" ,  /* Columns */
00059  "0028 0030" ,  /* Pixel spacing */
00060  "0028 0100" ,  /* Bits allocated */
00061  "0028 0101" ,  /* Bits stored */
00062  "0028 1052" ,  /* Rescale intercept */
00063  "0028 1053" ,  /* Rescale slope */
00064  "0028 1054" ,  /* Rescale type */
00065  "0028 0004" ,  /* Photometric interpretation */
00066  "0028 0103" ,  /* Pixel representation */
00067  "0028 0102" ,  /* High bit */
00068  "0028 1050" ,  /* Window center */
00069  "0028 1051" ,  /* Window width */
00070 
00071  "0008 0008" ,  /* ID Image type */
00072  "0008 0070" ,  /* ID Manufacturer */
00073  "0018 1310" ,  /* Acquisition Matrix */
00074 
00075  "0029 1010" ,  /* Siemens addendum #1 */
00076  "0029 1020" ,  /* Siemens addendum #2 */
00077 
00078 NULL } ;
00079 
00080 #define NUM_ELIST (sizeof(elist)/sizeof(char *)-1)
00081 
00082 #define E_SLICE_THICKNESS             0
00083 #define E_REPETITION_TIME             1
00084 #define E_SLICE_SPACING               2
00085 #define E_FIELD_OF_VIEW               3
00086 
00087 #define E_PATIENT_ORIENTATION         4
00088 #define E_IMAGE_POSITION              5
00089 #define E_IMAGE_ORIENTATION           6
00090 #define E_SLICE_LOCATION              7
00091 
00092 #define E_SAMPLES_PER_PIXEL           8
00093 #define E_NUMBER_OF_FRAMES            9
00094 #define E_ROWS                       10
00095 #define E_COLUMNS                    11
00096 #define E_PIXEL_SPACING              12
00097 #define E_BITS_ALLOCATED             13
00098 #define E_BITS_STORED                14
00099 #define E_RESCALE_INTERCEPT          15
00100 #define E_RESCALE_SLOPE              16
00101 #define E_RESCALE_TYPE               17
00102 #define E_PHOTOMETRIC_INTERPRETATION 18
00103 #define E_PIXEL_REPRESENTATION       19
00104 #define E_HIGH_BIT                   20
00105 #define E_WINDOW_CENTER              21
00106 #define E_WINDOW_WIDTH               22
00107 
00108 #define E_ID_IMAGE_TYPE              23    /* 28 Oct 2002: for Siemens mosaic */
00109 #define E_ID_MANUFACTURER            24
00110 #define E_ACQ_MATRIX                 25
00111 
00112 #define E_SIEMENS_1                  26    /* 31 Oct 2002 */
00113 #define E_SIEMENS_2                  27
00114 
00115 /*-----------------------------------------------------------------------------------*/
00116 /*! Read image(s) from a DICOM file, if possible.
00117 -------------------------------------------------------------------------------------*/
00118 
00119 MRI_IMARR * mri_read_dicom( char *fname )
00120 {
00121    char *ppp , *ddd ;
00122    off_t poff ;
00123    unsigned int plen ;
00124    char *epos[NUM_ELIST] ;
00125    int ii,jj , ee , bpp , datum ;
00126    int nx,ny,nz , swap , shift=0 ;
00127    float dx,dy,dz,dt ;
00128    MRI_IMARR *imar ;
00129    MRI_IMAGE *im ;
00130    void *iar ;
00131    FILE *fp ;
00132    short sbot,stop ;
00133    int have_orients=0 ;
00134    int ior,jor,kor ;
00135    static int nzoff=0 ;   /* for determining z-axis orientation/offset from multiple calls */
00136 
00137    int mosaic=0 , mos_nx,mos_ny , mos_ix,mos_iy,mos_nz ;  /* 28 Oct 2002 */
00138    Siemens_extra_info sexinfo ;                           /* 31 Oct 2002 */
00139    float xcen,ycen,zcen ;
00140    int use_xycen=0 ;
00141    float dxx,dyy,dzz ;
00142 
00143    char *eee ;
00144    float rescale_slope=0.0 , rescale_inter=0.0 ;  /* 23 Dec 2002 */
00145    float window_center=0.0 , window_width =0.0 ;
00146 
00147    char *sexi_start;   /* KRH 25 Jul 2003 */
00148    char *sexi_end;
00149 
00150 ENTRY("mri_read_dicom") ;
00151 
00152    if( str_sexinfo != NULL ){ free(str_sexinfo); str_sexinfo=NULL; }
00153 
00154    if( !mri_possibly_dicom(fname) ) RETURN(NULL) ;  /* 07 May 2003 */
00155 
00156    /* extract header info from file into a string
00157       - cf. mri_dicom_hdr.[ch]
00158       - run 'dicom_hdr -noname fname' to see the string format */
00159 
00160    mri_dicom_nohex(1) ;              /* don't print ints in hex mode */
00161    mri_dicom_noname(1) ;             /* don't print names, just tags */
00162    ppp = mri_dicom_header( fname ) ; /* print header to malloc()-ed string */
00163    if( ppp == NULL ) RETURN(NULL);   /* didn't work; not a DICOM file? */
00164 
00165    /* find out where the pixel array is in the file */
00166 
00167    mri_dicom_pxlarr( &poff , &plen ) ;
00168    if( plen <= 0 ){ free(ppp) ; RETURN(NULL); }
00169 
00170    /* check if file is actually this big (maybe it was truncated) */
00171 
00172    { unsigned int psiz , fsiz ;
00173      fsiz = THD_filesize( fname ) ;
00174      psiz = (unsigned int)(poff) + plen ;
00175      if( fsiz < psiz ){ free(ppp) ; RETURN(NULL); }
00176    }
00177 
00178    /* find positions in header of elements we care about */
00179 
00180    for( ee=0 ; ee < NUM_ELIST ; ee++ )
00181      epos[ee] = strstr(ppp,elist[ee]) ;
00182 
00183    /* see if the header has the elements we absolutely need */
00184 
00185    if( epos[E_ROWS]           == NULL ||
00186        epos[E_COLUMNS]        == NULL ||
00187        epos[E_BITS_ALLOCATED] == NULL   ){ free(ppp) ; RETURN(NULL); }
00188 
00189    /* check if we have 1 sample per pixel (can't deal with 3 or 4 now) */
00190 
00191    if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
00192       ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
00193       if( ddd == NULL ){ free(ppp) ; RETURN(NULL); }
00194       ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
00195       if( ii != 1 ){ free(ppp) ; RETURN(NULL); }
00196    }
00197 
00198    /* check if photometric interpretation is MONOCHROME (don't like PALETTE) */
00199 
00200    if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
00201       ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
00202       if( ddd == NULL ){ free(ppp) ; RETURN(NULL); }
00203    }
00204 
00205    /* check if we have 8, 16, or 32 bits per pixel */
00206 
00207    ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
00208    if( ddd == NULL ){ free(ppp); RETURN(NULL); }
00209    bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
00210    switch( bpp ){
00211       default: free(ppp); RETURN(NULL);    /* bad value */
00212       case  8: datum = MRI_byte ; break ;
00213       case 16: datum = MRI_short; break ;
00214       case 32: datum = MRI_int  ; break ;  /* probably not present in DICOM? */
00215    }
00216    bpp /= 8 ; /* now bytes per pixel, instead of bits */
00217 
00218    /*** Print some warnings if appropriate ***/
00219 
00220    /* check if BITS_STORED and HIGH_BIT are aligned */
00221 
00222 #define NWMAX 2
00223    if( epos[E_BITS_STORED] != NULL && epos[E_HIGH_BIT] != NULL ){
00224      int bs=0 , hb=0 ;
00225      ddd = strstr(epos[E_BITS_STORED],"//") ; sscanf(ddd+2,"%d",&bs) ;
00226      ddd = strstr(epos[E_HIGH_BIT],"//")    ; sscanf(ddd+2,"%d",&hb) ;
00227      if( bs != hb+1 ){
00228        static int nwarn=0 ;
00229        if( nwarn < NWMAX )
00230          fprintf(stderr,
00231                  "++ DICOM WARNING: file %s has Bits_Stored=%d and High_Bit=%d\n",
00232                  fname,bs,hb) ;
00233        if( nwarn == NWMAX )
00234          fprintf(stderr,"++ DICOM WARNING: no more Bits_Stored messages will be printed\n") ;
00235        nwarn++ ;
00236      }
00237    }
00238 
00239    /* check if Rescale is ordered */
00240    /* 23 Dec 2002: actually get the rescale params, if environment says to */
00241 
00242    eee = getenv("AFNI_DICOM_RESCALE") ;
00243    if( epos[E_RESCALE_INTERCEPT] != NULL && epos[E_RESCALE_SLOPE] != NULL ){
00244      if( eee == NULL || toupper(*eee) != 'Y' ){
00245        static int nwarn=0 ;
00246        if( nwarn < NWMAX )
00247          fprintf(stderr,
00248                  "++ DICOM WARNING: file %s has Rescale tags; setenv AFNI_DICOM_RESCALE YES to enforce them\n",
00249                  fname ) ;
00250        if( nwarn == NWMAX )
00251          fprintf(stderr,"++ DICOM WARNING: no more Rescale tags messages will be printed\n") ;
00252        nwarn++ ;
00253      } else {
00254        ddd = strstr(epos[E_RESCALE_INTERCEPT],"//") ; sscanf(ddd+2,"%f",&rescale_inter) ;
00255        ddd = strstr(epos[E_RESCALE_SLOPE    ],"//") ; sscanf(ddd+2,"%f",&rescale_slope) ;
00256      }
00257    }
00258 
00259    /* check if Window is ordered */
00260    /* 23 Dec 2002: actually get the window params, if environment says to */
00261 
00262    eee = getenv("AFNI_DICOM_WINDOW") ;
00263    if( epos[E_WINDOW_CENTER] != NULL && epos[E_WINDOW_WIDTH] != NULL ){
00264      if( eee == NULL || toupper(*eee) != 'Y' ){
00265        static int nwarn=0 ;
00266        if( nwarn < NWMAX )
00267          fprintf(stderr,
00268                  "++ DICOM WARNING: file %s has Window tags; setenv AFNI_DICOM_WINDOW YES to enforce them\n",
00269                  fname ) ;
00270        if( nwarn == NWMAX )
00271          fprintf(stderr,"++ DICOM WARNING: no more Window tags messages will be printed\n") ;
00272        nwarn++ ;
00273      } else {
00274        ddd = strstr(epos[E_WINDOW_CENTER],"//") ; sscanf(ddd+2,"%f",&window_center) ;
00275        ddd = strstr(epos[E_WINDOW_WIDTH ],"//") ; sscanf(ddd+2,"%f",&window_width ) ;
00276      }
00277    }
00278 
00279    /*** extract attributes of the image(s) to be read in ***/
00280 
00281    /* get image nx & ny */
00282 
00283    ddd = strstr(epos[E_ROWS],"//") ;                 /* 31 Oct 2002: */
00284    if( ddd == NULL ){ free(ppp) ; RETURN(NULL); }    /* Oops: ROWS is ny and */
00285    ny = 0 ; sscanf(ddd+2,"%d",&ny) ;                 /*       COLUMNS is nx! */
00286    if( ny < 2 ){ free(ppp) ; RETURN(NULL); }
00287 
00288    ddd = strstr(epos[E_COLUMNS],"//") ;
00289    if( ddd == NULL ){ free(ppp) ; RETURN(NULL); }
00290    nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
00291    if( nx < 2 ){ free(ppp) ; RETURN(NULL); }
00292 
00293    /* get number of slices */
00294 
00295    nz = 0 ;
00296    if( epos[E_NUMBER_OF_FRAMES] != NULL ){
00297      ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
00298      if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
00299    }
00300 
00301    /* if didn't get nz above, make up a value */
00302 
00303    if( nz == 0 ) nz = plen / (bpp*nx*ny) ;    /* compute from image array size */
00304    if( nz == 0 ){ free(ppp) ; RETURN(NULL); }
00305 
00306    /*-- 28 Oct 2002: Check if this is a Siemens mosaic.        --*/
00307    /*-- 02 Dec 2002: Don't use Acquisition Matrix anymore;
00308                      instead, use the Siemens extra info
00309                      in epos[E_SIEMENS_2].                     --*/
00310    /*-- 24 Dec 2002: Extract str_sexinfo even if not a mosaic. --*/
00311 
00312    if(        epos[E_ID_MANUFACTURER]            != NULL &&
00313        strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
00314               epos[E_SIEMENS_2]                  != NULL    ){
00315 
00316      int len=0,loc=0 , aa,bb ;
00317      sscanf(epos[E_SIEMENS_2],"%x%x%d [%d" , &aa,&bb , &len,&loc ) ;
00318      if( len > 0 && loc > 0 ){
00319        fp = fopen( fname , "rb" ) ;
00320        if( fp != NULL ){
00321          str_sexinfo = extract_bytes_from_file( fp, (off_t)loc, (size_t)len, 1 ) ;
00322          fclose(fp) ;
00323        }
00324      }
00325    }
00326 
00327    /*-- process str_sexinfo only if this is marked as a mosaic image --*/
00328 
00329    if(        epos[E_ID_IMAGE_TYPE]              != NULL &&
00330        strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC")    != NULL &&
00331        str_sexinfo                               != NULL   ){
00332 
00333      /* 31 Oct 2002: extract extra Siemens info from str_sexinfo */
00334 
00335      /* KRH 25 Jul 2003 if start and end markers are present for
00336       * Siemens extra info, cut string down to those boundaries */
00337 
00338      sexi_start = strstr(str_sexinfo, "### ASCCONV BEGIN ###");
00339      sexi_end = strstr(str_sexinfo, "### ASCCONV END ###");
00340      if ((sexi_start != NULL) && (sexi_end != NULL)) {
00341        char *sexi_tmp;
00342        int sexi_size;
00343 
00344        sexi_size = sexi_end - sexi_start + 19 ;
00345        sexi_tmp = AFMALL( char, sexi_size );
00346        memcpy(sexi_tmp,sexi_start,sexi_size);
00347        free(str_sexinfo);
00348        str_sexinfo = sexi_tmp;
00349      }
00350      /* end KRH 25 Jul 2003 change */
00351 
00352      sexinfo.good = 0 ;  /* start by marking it as bad */
00353      for(ii = 0; ii < 3; ii++) sexinfo.have_data[ii] = 0; /* 25 Feb 03 Initialize new member KRH */
00354 
00355      get_siemens_extra_info( str_sexinfo , &sexinfo ) ;
00356 
00357      if( sexinfo.good ){                                   /* if data is good */
00358 
00359        /* compute size of mosaic layout
00360           as 1st integer whose square is >= # of images in mosaic */
00361 
00362        for( mos_ix=1 ; mos_ix*mos_ix < sexinfo.mosaic_num ; mos_ix++ ) ; /* nada */
00363        mos_iy = mos_ix ;
00364 
00365        mos_nx = nx / mos_ix ; mos_ny = ny / mos_iy ;  /* sub-image dimensions */
00366 
00367        if( mos_ix*mos_nx == nx &&               /* Sub-images must fit nicely */
00368            mos_iy*mos_ny == ny    ){            /* into super-image layout.   */
00369 
00370            static int nwarn=0 ;
00371 
00372            /* should be tagged as a 1 slice image thus far */
00373 
00374            if( nz > 1 ){
00375              fprintf(stderr,
00376                      "** DICOM ERROR: %dx%d Mosaic of %dx%d images in file %s, but also have nz=%d\n",
00377                      mos_ix,mos_iy,mos_nx,mos_ny,fname,nz) ;
00378              free(ppp) ; RETURN(NULL) ;
00379            }
00380 
00381            /* mark as a mosaic */
00382 
00383            mosaic = 1 ;
00384            mos_nz = mos_ix * mos_iy ;   /* number of slices in mosaic */
00385            if( nwarn < NWMAX )
00386              fprintf(stderr,"++ DICOM NOTICE: %dx%d Siemens Mosaic of %d %dx%d images in file %s\n",
00387                     mos_ix,mos_iy,sexinfo.mosaic_num,mos_nx,mos_ny,fname) ;
00388            if( nwarn == NWMAX )
00389              fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosiac messages will be printed\n") ;
00390            nwarn++ ;
00391 
00392        } /* end of if mosaic sizes are reasonable */
00393 
00394        else {                        /* warn about bad mosaic sizes */
00395          static int nwarn=0 ;
00396          if( nwarn < NWMAX )
00397            fprintf(stderr,
00398                    "++ DICOM WARNING: bad Siemens Mosaic params: nx=%d ny=%d ix=%d iy=%d imx=%d imy=%d\n",
00399                    mos_nx,mos_ny , mos_ix,mos_iy , nx,ny ) ;
00400          if( nwarn == NWMAX )
00401            fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic param messages will be printed\n");
00402          nwarn++ ;
00403        }
00404 
00405      } /* end of if sexinfo was good */
00406 
00407      else {                  /* warn if sexinfo was bad */
00408        static int nwarn=0 ;
00409        if( nwarn < NWMAX )
00410          fprintf(stderr,"++ DICOM WARNING: indecipherable Siemens Mosaic info (%s) in file %s\n",
00411                  elist[E_SIEMENS_2] , fname ) ;
00412        if( nwarn == NWMAX )
00413          fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic info messages will be printed\n");
00414        nwarn++ ;
00415      }
00416 
00417    } /* end of if a Siemens mosaic */
00418 
00419    /*-- try to get dx, dy, dz, dt --*/
00420 
00421    dx = dy = dz = dt = 0.0 ;
00422 
00423    /* dx,dy first */
00424 
00425    if( epos[E_PIXEL_SPACING] != NULL ){
00426      ddd = strstr(epos[E_PIXEL_SPACING],"//") ;
00427      if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00428      if( dy == 0.0 && dx > 0.0 ) dy = dx ;
00429    }
00430    if( dx == 0.0 && epos[E_FIELD_OF_VIEW] != NULL ){
00431      ddd = strstr(epos[E_FIELD_OF_VIEW],"//") ;
00432      if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00433      if( dx > 0.0 ){
00434        if( dy == 0.0 ) dy = dx ;
00435        dx /= nx ; dy /= ny ;
00436      }
00437    }
00438 
00439    /*-- 27 Nov 2002: fix stupid GE error,
00440                      where the slice spacing is really the slice gap --*/
00441 
00442    { int stupid_ge_fix , no_stupidity ;
00443      float sp=0.0 , th=0.0 ;
00444      static int nwarn=0 ;
00445 
00446      eee           = getenv("AFNI_SLICE_SPACING_IS_GAP") ;
00447      stupid_ge_fix = (eee != NULL && (*eee=='Y' || *eee=='y') ) ;
00448      no_stupidity  = (eee != NULL && (*eee=='N' || *eee=='n') ) ;  /* 03 Mar 2003 */
00449 
00450      if( epos[E_SLICE_SPACING] != NULL ){                  /* get reported slice spacing */
00451        ddd = strstr(epos[E_SLICE_SPACING],"//") ;
00452        if( ddd != NULL ) sscanf( ddd+2 , "%f" , &sp ) ;
00453      }
00454      if( epos[E_SLICE_THICKNESS] != NULL ){                /* get reported slice thickness */
00455        ddd = strstr(epos[E_SLICE_THICKNESS],"//") ;
00456        if( ddd != NULL ) sscanf( ddd+2 , "%f" , &th ) ;
00457      }
00458 
00459      th = fabs(th) ; sp = fabs(sp) ;                       /* we don't use the sign */
00460 
00461      if( stupid_ge_fix ){                                  /* always be stupid */
00462        dz = sp+th ;
00463      } else {
00464 
00465        if( no_stupidity && sp > 0.0 )                      /* 13 Jan 2004: if 'NO', then */
00466          dz = sp ;                                         /* always use spacing if present */
00467        else
00468          dz = (sp > th) ? sp : th ;                        /* the correct-ish DICOM way */
00469 
00470 #define GFAC 0.99
00471 
00472        if( !no_stupidity ){                                /* unless stupidity is turned off */
00473          if( sp > 0.0 && sp < GFAC*th ) dz = sp+th ;       /* the stupid GE way again */
00474 
00475          if( sp > 0.0 && sp < GFAC*th && nwarn < NWMAX ){
00476            fprintf(stderr,
00477                    "++ DICOM WARNING: file %s has Slice_Spacing=%f smaller than Slice_Thickness=%f\n",
00478                    fname , sp , th ) ;
00479            if( nwarn == 0 )
00480             fprintf(stderr,
00481               "\n"
00482               "++  Setting environment variable AFNI_SLICE_SPACING_IS_GAP       ++\n"
00483               "++   to YES will make the center-to-center slice distance        ++\n"
00484               "++   be set to Slice_Spacing+Slice_Thickness=%6.3f.             ++\n"
00485               "++  This is against the DICOM standard [attribute (0018,0088)    ++\n"
00486               "++   is defined as the center-to-center spacing between slices,  ++\n"
00487               "++   NOT as the edge-to-edge gap between slices], but it seems   ++\n"
00488               "++   to be necessary for some GE scanners.                       ++\n"
00489               "++                                                               ++\n"
00490               "++  This correction has been made on this data: dz=%6.3f.       ++\n"
00491               "++                                                               ++\n"
00492               "++  Setting AFNI_SLICE_SPACING_IS_GAP to NO means that the       ++\n"
00493               "++  DICOM Slice_Spacing variable will be used for dz, replacing  ++\n"
00494               "++  the Slice_Thickness variable.  This usage may be required    ++\n"
00495               "++  for some pulse sequences on Phillips scanners.               ++\n"
00496               "\n\a" ,
00497              sp+th , dz ) ;
00498          }
00499          if( sp > 0.0 && sp < th && nwarn == NWMAX )
00500            fprintf(stderr,"++ DICOM WARNING: no more Slice_Spacing messages will be printed\n") ;
00501          nwarn++ ;
00502        }
00503      }
00504      if( dz == 0.0 && dx != 0.0 ) dz = 1.0 ;               /* nominal dz */
00505 
00506    } /*-- end of dz code, with all its stupidities --*/
00507 
00508    /* get dt */
00509 
00510    if( epos[E_REPETITION_TIME] != NULL ){
00511      ddd = strstr(epos[E_REPETITION_TIME],"//") ;
00512      if( ddd != NULL ){
00513        sscanf( ddd+2 , "%f" , &dt ) ;
00514        dt *= 0.001 ;   /* ms to s */
00515      }
00516    }
00517 
00518    /* check if we might have 16 bit unsigned data that fills all bits */
00519 
00520 #if 0
00521    if( bpp == 2 ){
00522      if( epos[E_PIXEL_REPRESENTATION] != NULL ){
00523        ddd = strstr(epos[E_PIXEL_REPRESENTATION],"//") ;
00524        if( ddd != NULL ){
00525          ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
00526          if( ii == 1 ) shift = -1 ;
00527        }
00528      }
00529      if( shift == 0 && epos[E_HIGH_BIT] != NULL ){
00530        ddd = strstr(epos[E_HIGH_BIT],"//") ;
00531        if( ddd != NULL ){
00532          ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
00533          if( ii == 15 ) shift = 1 ;
00534        }
00535      }
00536      sbot = 32767 ; stop = -32767 ;
00537    }
00538 #endif
00539 
00540    /** Finally! Read images from file. **/
00541 
00542    fp = fopen( fname , "rb" ) ;
00543    if( fp == NULL ){ free(ppp) ; RETURN(NULL); }
00544    lseek( fileno(fp) , poff , SEEK_SET ) ;
00545 
00546    INIT_IMARR(imar) ;
00547 
00548    /* DICOM files are stored in LSB first (little endian) mode */
00549 
00550    RWC_set_endianosity() ; swap = !LITTLE_ENDIAN_ARCHITECTURE ;
00551 
00552    /* 28 Oct 2002: must allow for 2D mosaic mode */
00553 
00554    if( !mosaic ){   /*-- 28 Oct 2002: old method, not a mosaic --*/
00555 
00556     for( ii=0 ; ii < nz ; ii++ ){
00557       im = mri_new( nx , ny , datum ) ;    /* new MRI_IMAGE struct */
00558       iar = mri_data_pointer( im ) ;       /* data array in struct */
00559       fread( iar , bpp , nx*ny , fp ) ;    /* read data directly into it */
00560 
00561       if( swap ){                          /* swap bytes? */
00562         switch( im->pixel_size ){
00563           default: break ;
00564           case 2:  swap_twobytes (   im->nvox, iar ) ; break ;  /* short */
00565           case 4:  swap_fourbytes(   im->nvox, iar ) ; break ;  /* int, float */
00566           case 8:  swap_fourbytes( 2*im->nvox, iar ) ; break ;  /* complex */
00567         }
00568         im->was_swapped = 1 ;
00569       }
00570 
00571 #if 0
00572       if( shift == 1 ){
00573         switch( datum ){
00574           case MRI_short:{
00575             short * sar = (short *) iar ;
00576             for( jj=0 ; jj < im->nvox ; jj++ ){
00577               sbot = MIN( sar[jj] , sbot ) ;
00578               stop = MAX( sar[jj] , stop ) ;
00579             }
00580           }
00581           break ;
00582         }
00583       }
00584 #endif
00585 
00586       /* store auxiliary data in image struct */
00587 
00588       if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
00589         im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
00590       }
00591       if( dt > 0.0 ) im->dt = dt ;
00592 
00593       ADDTO_IMARR(imar,im) ;
00594     }
00595 
00596    } else {   /*-- 28 Oct 2002:  is a 2D mosaic --*/
00597 
00598      char *dar , *iar ;
00599      int last_ii=-1 , nvox , yy,xx,nxx ;
00600 
00601      nvox = mos_nx*mos_ny*mos_nz ;         /* total number of voxels */
00602      dar  = (char*)calloc(bpp,nvox) ;            /* make space for super-image */
00603      fread( dar , bpp , nvox , fp ) ;    /* read data directly into it */
00604      if( swap ){                        /* swap bytes? */
00605        switch( bpp ){
00606          default: break ;
00607          case 2:  swap_twobytes (   nvox, dar ) ; break ;  /* short */
00608          case 4:  swap_fourbytes(   nvox, dar ) ; break ;  /* int, float */
00609          case 8:  swap_fourbytes( 2*nvox, dar ) ; break ;  /* complex */
00610        }
00611      }
00612 
00613      /* load data from dar into images */
00614 
00615      nxx = mos_nx * mos_ix ;              /* # pixels per mosaic line */
00616 
00617      for( yy=0 ; yy < mos_iy ; yy++ ){    /* loop over sub-images in mosaic */
00618        for( xx=0 ; xx < mos_ix ; xx++ ){
00619 
00620          im = mri_new( mos_nx , mos_ny , datum ) ;
00621          iar = mri_data_pointer( im ) ;             /* sub-image array */
00622 
00623          /* copy data rows from dar into iar */
00624 
00625          for( jj=0 ; jj < mos_ny ; jj++ )  /* loop over rows inside sub-image */
00626            memcpy( iar + jj*mos_nx*bpp ,
00627                    dar + xx*mos_nx*bpp + (jj+yy*mos_ny)*nxx*bpp ,
00628                    mos_nx*bpp                                    ) ;
00629 
00630          if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
00631            im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
00632          }
00633          if( dt > 0.0 ) im->dt = dt ;
00634          if( swap ) im->was_swapped = 1 ;
00635 
00636          ADDTO_IMARR(imar,im) ;
00637        }
00638      }
00639      free(dar) ;  /* don't need no more; copied all data out of it now */
00640 
00641      /* truncate zero images out of tail of mosaic */
00642 
00643      if( sexinfo.good ){  /* the new way: use the mosaic count from Siemens extra info */
00644 
00645        if( sexinfo.mosaic_num < IMARR_COUNT(imar) )
00646          TRUNCATE_IMARR(imar,sexinfo.mosaic_num) ;
00647 
00648      } else {  /* the old way: find the last image with a nonzero value inside */
00649 
00650        for( ii=mos_nz-1 ; ii >= 0 ; ii-- ){  /* loop backwards over images */
00651          im = IMARR_SUBIM(imar,ii) ;
00652          switch( im->kind ){
00653            case MRI_short:{
00654              short *ar = mri_data_pointer( im ) ;
00655              for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
00656              if( jj < im->nvox ) last_ii = ii ;
00657            }
00658            break ;
00659 
00660            case MRI_byte:{
00661              byte *ar = mri_data_pointer( im ) ;
00662              for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
00663              if( jj < im->nvox ) last_ii = ii ;
00664            }
00665            break ;
00666 
00667            case MRI_int:{
00668              int *ar = mri_data_pointer( im ) ;
00669              for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
00670              if( jj < im->nvox ) last_ii = ii ;
00671            }
00672            break ;
00673          }
00674          if( last_ii >= 0 ) break ;
00675        }
00676 
00677        if( last_ii <= 0 ) last_ii = 1 ;
00678        if( last_ii+1 < IMARR_COUNT(imar) ) TRUNCATE_IMARR(imar,last_ii+1) ;
00679 
00680      } /* end of truncating off all zero images at end */
00681 
00682 #if 0
00683 fprintf(stderr,"\nmri_read_dicom Mosaic: mos_nx=%d mos_ny=%d mos_ix=%d mos_iy=%d slices=%d\n",
00684 mos_nx,mos_ny,mos_ix,mos_iy,IMARR_COUNT(imar)) ;
00685 MCHECK ;
00686 #endif
00687 
00688    } /* end of mosaic input mode */
00689 
00690    fclose(fp) ;     /* 10 Sep 2002: oopsie - forgot to close file */
00691 
00692    /*-- 23 Dec 2002: implement Rescale, if ordered --*/
00693 
00694    if( rescale_slope > 0.0 ){
00695      for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
00696        im = IMARR_SUBIM(imar,ii) ;
00697        switch( im->kind ){
00698          case MRI_byte:{
00699            byte *ar = mri_data_pointer( im ) ;
00700            for( jj=0 ; jj < im->nvox ; jj++ )
00701              ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00702          }
00703          break ;
00704 
00705          case MRI_short:{
00706            short *ar = mri_data_pointer( im ) ;
00707            for( jj=0 ; jj < im->nvox ; jj++ )
00708              ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00709          }
00710          break ;
00711 
00712          case MRI_int:{
00713            int *ar = mri_data_pointer( im ) ;
00714            for( jj=0 ; jj < im->nvox ; jj++ )
00715              ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00716          }
00717          break ;
00718        }
00719      }
00720    } /* end of Rescale */
00721 
00722    /*-- 23 Dec 2002: implement Window, if ordered --*/
00723    /*                section C.11.2.1.2 (page 503)  */
00724 
00725    if( window_width >= 1.0 ){
00726      float wbot,wtop,wfac ;
00727      int ymax=0 ;
00728 
00729      /* get output range */
00730 
00731      ddd = strstr(epos[E_BITS_STORED],"//") ;
00732      if( ddd != NULL ){
00733        ymax = 0 ; sscanf(ddd+2,"%d",&ymax) ;
00734        if( ymax > 0 ) ymax = (1 << ymax) - 1 ;
00735      }
00736      if( ymax <= 0 ){
00737        switch( IMARR_SUBIM(imar,0)->kind ){
00738          case MRI_byte:  ymax = MRI_maxbyte  ; break ;
00739          case MRI_short: ymax = MRI_maxshort ; break ;
00740          case MRI_int:   ymax = MRI_maxint   ; break ;
00741        }
00742      }
00743                                           /** window_width == 1 is special **/
00744      if( window_width == 1.0 ){           /** binary threshold case **/
00745 
00746        wbot = window_center - 0.5 ;       /* the threshold */
00747 
00748        for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
00749          im = IMARR_SUBIM(imar,ii) ;
00750          switch( im->kind ){
00751            case MRI_byte:{
00752              byte *ar = mri_data_pointer( im ) ;
00753              for( jj=0 ; jj < im->nvox ; jj++ )
00754                ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
00755            }
00756            break ;
00757   
00758            case MRI_short:{
00759              short *ar = mri_data_pointer( im ) ;
00760              for( jj=0 ; jj < im->nvox ; jj++ )
00761                ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
00762            }
00763            break ;
00764   
00765            case MRI_int:{
00766              int *ar = mri_data_pointer( im ) ;
00767              for( jj=0 ; jj < im->nvox ; jj++ )
00768                ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
00769            }
00770            break ;
00771          }
00772        }
00773 
00774      } else {                             /** linear windowing case **/
00775 
00776        wbot = (window_center - 0.5) - 0.5*(window_width-1.0) ;
00777        wtop = (window_center - 0.5) + 0.5*(window_width-1.0) ;
00778        wfac = ymax                  /     (window_width-1.0) ;
00779 
00780        for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
00781          im = IMARR_SUBIM(imar,ii) ;
00782          switch( im->kind ){
00783            case MRI_byte:{
00784              byte *ar = mri_data_pointer( im ) ;
00785              for( jj=0 ; jj < im->nvox ; jj++ )
00786                ar[jj] = (ar[jj] <= wbot) ? 0
00787                        :(ar[jj] >  wtop) ? ymax
00788                                          : wfac*(ar[jj]-wbot)+0.499 ;
00789            }
00790            break ;
00791 
00792            case MRI_short:{
00793              short *ar = mri_data_pointer( im ) ;
00794              for( jj=0 ; jj < im->nvox ; jj++ )
00795                ar[jj] = (ar[jj] <= wbot) ? 0
00796                        :(ar[jj] >  wtop) ? ymax
00797                                          : wfac*(ar[jj]-wbot)+0.499 ;
00798            }
00799            break ;
00800 
00801            case MRI_int:{
00802              int *ar = mri_data_pointer( im ) ;
00803              for( jj=0 ; jj < im->nvox ; jj++ )
00804                ar[jj] = (ar[jj] <= wbot) ? 0
00805                        :(ar[jj] >  wtop) ? ymax
00806                                          : wfac*(ar[jj]-wbot)+0.499 ;
00807            }
00808            break ;
00809          }
00810        }
00811      }
00812    } /* end of Window */
00813 
00814    /*-- store some extra information in MRILIB globals, too? --*/
00815 
00816    if( dt > 0.0 && MRILIB_tr <= 0.0 ) MRILIB_tr = dt ;  /* TR */
00817 
00818    /*-- try to get image orientation fields (also, set ior,jor,kor) --*/
00819 
00820    if( epos[E_IMAGE_ORIENTATION] != NULL ){    /* direction cosines of image plane */
00821 
00822      ddd = strstr(epos[E_IMAGE_ORIENTATION],"//") ;
00823      if( ddd != NULL ){
00824        float xc1=0.0,xc2=0.0,xc3=0.0 , yc1=0.0,yc2=0.0,yc3=0.0 ;
00825        float xn,yn ; int qq ;
00826        qq = sscanf(ddd+2,"%f\\%f\\%f\\%f\\%f\\%f",&xc1,&xc2,&xc3,&yc1,&yc2,&yc3) ;
00827        xn = sqrt( xc1*xc1 + xc2*xc2 + xc3*xc3 ) ; /* vector norms */
00828        yn = sqrt( yc1*yc1 + yc2*yc2 + yc3*yc3 ) ;
00829        if( qq == 6 && xn > 0.0 && yn > 0.0 ){     /* both vectors OK */
00830 
00831          xc1 /= xn ; xc2 /= xn ; xc3 /= xn ;      /* normalize vectors */
00832          yc1 /= yn ; yc2 /= yn ; yc3 /= yn ;
00833 
00834          if( !use_MRILIB_xcos ){
00835            MRILIB_xcos[0] = xc1 ; MRILIB_xcos[1] = xc2 ;  /* save direction */
00836            MRILIB_xcos[2] = xc3 ; use_MRILIB_xcos = 1 ;   /* cosine vectors */
00837          }
00838 
00839          if( !use_MRILIB_ycos ){
00840            MRILIB_ycos[0] = yc1 ; MRILIB_ycos[1] = yc2 ;
00841            MRILIB_ycos[2] = yc3 ; use_MRILIB_ycos = 1 ;
00842          }
00843 
00844          /* x-axis orientation */
00845          /* ior determines which spatial direction is x-axis  */
00846          /* and is the direction that has the biggest change */
00847 
00848          dxx = fabs(xc1) ; ior = 1 ;
00849          dyy = fabs(xc2) ; if( dyy > dxx ){ ior=2; dxx=dyy; }
00850          dzz = fabs(xc3) ; if( dzz > dxx ){ ior=3;        }
00851          dxx = MRILIB_xcos[ior-1] ; if( dxx < 0. ) ior = -ior;
00852 
00853          if( MRILIB_orients[0] == '\0' ){
00854            switch( ior ){
00855              case -1: MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; break;
00856              case  1: MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; break;
00857              case -2: MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; break;
00858              case  2: MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; break;
00859              case  3: MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; break;
00860              case -3: MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; break;
00861              default: MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; break;
00862            }
00863          }
00864 
00865          /* y-axis orientation */
00866          /* jor determines which spatial direction is y-axis  */
00867          /* and is the direction that has the biggest change */
00868 
00869          dxx = fabs(yc1) ; jor = 1 ;
00870          dyy = fabs(yc2) ; if( dyy > dxx ){ jor=2; dxx=dyy; }
00871          dzz = fabs(yc3) ; if( dzz > dxx ){ jor=3;        }
00872          dyy = MRILIB_ycos[jor-1] ; if( dyy < 0. ) jor = -jor;
00873          if( MRILIB_orients[2] == '\0' ){
00874            switch( jor ){
00875              case -1: MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; break;
00876              case  1: MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; break;
00877              case -2: MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; break;
00878              case  2: MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; break;
00879              case  3: MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; break;
00880              case -3: MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; break;
00881              default: MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; break;
00882            }
00883          }
00884 
00885          MRILIB_orients[6] = '\0' ;   /* terminate orientation string */
00886 
00887          kor = 6 - abs(ior)-abs(jor) ;   /* which spatial direction is z-axis */
00888                                          /* where 1=L-R, 2=P-A, 3=I-S */
00889          have_orients = 1 ;
00890 #if 0
00891 fprintf(stderr,"MRILIB_orients=%s (from IMAGE_ORIENTATION)\n",MRILIB_orients) ;
00892 #endif
00893        }
00894      }
00895 
00896    } else if( epos[E_PATIENT_ORIENTATION] != NULL ){  /* symbolic orientation of image */
00897                                                       /* [not so useful, or common] */
00898      ddd = strstr(epos[E_PATIENT_ORIENTATION],"//") ;
00899      if( ddd != NULL ){
00900        char xc='\0' , yc='\0' ;
00901        sscanf(ddd+2,"%c\\%c",&xc,&yc) ;   /* e.g., "L\P" */
00902        switch( toupper(xc) ){
00903          case 'L': MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; ior=-1; break;
00904          case 'R': MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; ior= 1; break;
00905          case 'P': MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; ior=-2; break;
00906          case 'A': MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; ior= 2; break;
00907          case 'F': MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; ior= 3; break;  /* F = foot */
00908          case 'H': MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; ior=-3; break;  /* H = head */
00909          default:  MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; ior= 0; break;
00910        }
00911        switch( toupper(yc) ){
00912          case 'L': MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; jor=-1; break;
00913          case 'R': MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; jor= 1; break;
00914          case 'P': MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; jor=-2; break;
00915          case 'A': MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; jor= 2; break;
00916          case 'F': MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; jor= 3; break;
00917          case 'H': MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; jor=-3; break;
00918          default:  MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; jor= 0; break;
00919        }
00920        MRILIB_orients[6] = '\0' ;      /* terminate orientation string */
00921        kor = 6 - abs(ior)-abs(jor) ;   /* which spatial direction is z-axis */
00922        have_orients = (ior != 0 && jor != 0) ;
00923      }
00924 
00925    }  /* end of 2D image orientation */
00926 
00927    /*-- try to get image offset (position), if have orientation from above --*/
00928 
00929    if( nzoff == 0 && have_orients && mosaic && sexinfo.good ){  /* 01 Nov 2002: use Siemens mosaic info */
00930      int qq ;
00931      float z0, z1 ;
00932      /* 25 Feb 2003 changing error checking for mosaics missing one or more *
00933       * dimension of slice coordinates                                 KRH  */
00934      if (sexinfo.have_data[kor-1]) {
00935        z0 = sexinfo.slice_xyz[0][kor-1] ;   /* kor from orients above */
00936        z1 = sexinfo.slice_xyz[1][kor-1] ;   /* z offsets of 1st 2 slices */
00937      } else {                  /* warn if sexinfo was bad */
00938        static int nwarn=0 ;
00939        if( nwarn < NWMAX )
00940          fprintf(stderr,"++ DICOM WARNING: Unusable coord. in Siemens Mosaic info (%s) in file %s\n",
00941                  elist[E_SIEMENS_2] , fname ) ;
00942        if( nwarn == NWMAX )
00943          fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic info messages will be printed\n");
00944        nwarn++ ;
00945      }
00946 
00947 
00948 #if 0
00949      /* Save x,y center of this 1st slice */
00950 
00951      xcen = sexinfo.slice_xyz[0][abs(ior)-1] ;
00952      ycen = sexinfo.slice_xyz[0][abs(jor)-1] ; use_xycen = 1 ;
00953 #endif
00954 
00955      /* finish z orientation now */
00956 
00957      if( z1-z0 < 0.0 ) kor = -kor ;     /* reversed orientation */
00958      if( MRILIB_orients[4] == '\0' ){
00959        switch( kor ){
00960          case  1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
00961          case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
00962          case  2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
00963          case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
00964          case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
00965          case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
00966          default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
00967        }
00968      }
00969      MRILIB_orients[6] = '\0' ;
00970 
00971 #if 0
00972 fprintf(stderr,"z0=%f z1=%f kor=%d MRILIB_orients=%s\n",z0,z1,kor,MRILIB_orients) ;
00973 #endif
00974 
00975      MRILIB_zoff = z0 ; use_MRILIB_zoff = 1 ;
00976      if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
00977    }
00978 
00979    /** use image position vector to set offsets,
00980        and (2cd time in) the z-axis orientation **/
00981 
00982    if( nzoff < 2 && epos[E_IMAGE_POSITION] != NULL && have_orients ){
00983      ddd = strstr(epos[E_IMAGE_POSITION],"//") ;
00984      if( ddd != NULL ){
00985        float xyz[3] ; int qq ;
00986        qq = sscanf(ddd+2,"%f\\%f\\%f",xyz,xyz+1,xyz+2) ;
00987        if( qq == 3 ){
00988          static float zoff ;      /* saved from nzoff=0 case */
00989          float zz = xyz[kor-1] ;  /* kor from orients above */
00990 
00991 #if 0
00992 fprintf(stderr,"IMAGE_POSITION=%f %f %f  kor=%d\n",xyz[0],xyz[1],xyz[2],kor) ;
00993 #endif
00994 
00995          if( nzoff == 0 ){  /* 1st DICOM image */
00996 
00997            zoff = zz ;      /* save this for 2nd image calculation */
00998 
00999            /* 01 Nov 2002: in mosaic case, may have set this already */
01000 
01001            if( MRILIB_orients[4] == '\0' ){
01002              switch( kor ){   /* may be changed on second image */
01003                case  1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
01004                case  2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
01005                case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
01006                default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
01007              }
01008              MRILIB_orients[6] = '\0' ;
01009            }
01010 
01011            /* Save x,y offsets of this 1st slice */
01012 
01013            qq = abs(ior) ;
01014            MRILIB_xoff = xyz[qq-1] ; use_MRILIB_xoff = 1 ;
01015            if( ior > 0 ) MRILIB_xoff = -MRILIB_xoff ;
01016 
01017            qq = abs(jor) ;
01018            MRILIB_yoff = xyz[qq-1] ; use_MRILIB_yoff = 1 ;
01019            if( jor > 0 ) MRILIB_yoff = -MRILIB_yoff ;
01020 
01021            /* 01 Nov 2002: adjust x,y offsets for mosaic */
01022 
01023            if( mosaic ){
01024              if( MRILIB_xoff < 0.0 ) MRILIB_xoff += 0.5*dx*mos_nx*(mos_ix-1) ;
01025              else                    MRILIB_xoff -= 0.5*dx*mos_nx*(mos_ix-1) ;
01026              if( MRILIB_yoff < 0.0 ) MRILIB_yoff += 0.5*dy*mos_ny*(mos_iy-1) ;
01027              else                    MRILIB_yoff -= 0.5*dy*mos_ny*(mos_iy-1) ;
01028            }
01029 
01030          } else if( nzoff == 1 && !use_MRILIB_zoff ){  /* 2nd DICOM image */
01031 
01032            float qoff = zz - zoff ;    /* vive la difference */
01033            if( qoff < 0 ) kor = -kor ; /* kor determines z-axis orientation */
01034 #if 0
01035 fprintf(stderr,"  nzoff=1 kor=%d qoff=%f\n",kor,qoff) ;
01036 #endif
01037            switch( kor ){
01038              case  1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
01039              case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
01040              case  2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
01041              case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
01042              case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
01043              case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
01044              default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
01045            }
01046            MRILIB_orients[6] = '\0' ;
01047 
01048            /* save spatial offset of first slice              */
01049            /* [this needs to be positive in the direction of] */
01050            /* [the -z axis, so may need to change its sign  ] */
01051 
01052            MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
01053            if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
01054          }
01055          nzoff++ ;  /* 3rd and later images don't count for z-orientation */
01056        }
01057      }
01058    }  /* end of using image position */
01059 
01060    /** 23 Dec 2002:
01061        use slice location value to set z-offset,
01062        and (2cd time in) the z-axis orientation
01063        -- only try this if image position vector (code above) isn't present
01064           AND if we don't have a mosaic image (which already did this stuff)
01065        -- shouldn't be necessary, since slice location is deprecated        **/
01066 
01067    else if( nzoff < 2 && epos[E_SLICE_LOCATION] != NULL && have_orients && !mosaic ){
01068      ddd = strstr(epos[E_SLICE_LOCATION],"//") ;
01069      if( ddd != NULL ){
01070        float zz ; int qq ;
01071        qq = sscanf(ddd+2,"%f",&zz) ;
01072        if( qq == 1 ){
01073          static float zoff ;      /* saved from nzoff=0 case */
01074 
01075 #if 0
01076 fprintf(stderr,"SLICE_LOCATION = %f\n",zz) ;
01077 #endif
01078 
01079          if( nzoff == 0 ){  /* 1st DICOM image */
01080 
01081            zoff = zz ;      /* save this for 2nd image calculation */
01082 
01083            if( MRILIB_orients[4] == '\0' ){
01084              switch( kor ){   /* may be changed on second image */
01085                case  1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
01086                case  2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
01087                case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
01088                default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
01089              }
01090              MRILIB_orients[6] = '\0' ;
01091            }
01092 
01093          } else if( nzoff == 1 && !use_MRILIB_zoff ){  /* 2nd DICOM image */
01094 
01095            float qoff = zz - zoff ;    /* vive la difference */
01096            if( qoff < 0 ) kor = -kor ; /* kor determines z-axis orientation */
01097            switch( kor ){
01098              case  1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
01099              case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
01100              case  2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
01101              case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
01102              case  3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
01103              case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
01104              default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
01105            }
01106            MRILIB_orients[6] = '\0' ;
01107 
01108            /* save spatial offset of first slice              */
01109            /* [this needs to be positive in the direction of] */
01110            /* [the -z axis, so may need to change its sign  ] */
01111 
01112            MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
01113            if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
01114          }
01115          nzoff++ ;  /* 3rd and later images don't count for z-orientation */
01116        }
01117      }
01118    } /* end of using slice location */
01119 
01120    /* perhaps shift data shorts */
01121 
01122 #if 0
01123    if( shift == 1 ){
01124      switch( datum ){
01125        case MRI_short:{
01126          if( sbot < 0 ){
01127            unsigned short *sar ; int nvox = IMARR_SUBIM(imar,0)->nvox ;
01128            for( ii=0 ; ii < nz ; ii++ ){
01129              sar = mri_data_pointer( IMARR_SUBIM(imar,ii) ) ;
01130              for( jj=0 ; jj < nvox ; jj++ ) sar[jj] >>= 1 ;
01131            }
01132          }
01133        }
01134        break ;
01135      }
01136    }
01137 #endif
01138 
01139    RETURN( imar );
01140 }
01141 
01142 /*------------------------------------------------------------------------------*/
01143 /*! Count images in a DICOM file, if possible.
01144 --------------------------------------------------------------------------------*/
01145 
01146 int mri_imcount_dicom( char *fname )
01147 {
01148    char *ppp , *ddd ;
01149    off_t poff ;
01150    unsigned int plen ;
01151    char *epos[NUM_ELIST] ;
01152    int ii , ee , bpp , datum ;
01153    int nx,ny,nz ;
01154 
01155    int mosaic=0 , mos_nx,mos_ny , mos_ix,mos_iy,mos_nz ;  /* 28 Oct 2002 */
01156    Siemens_extra_info sexinfo ;                           /* 02 Dec 2002 */
01157    
01158    char *sexi_start;   /* KRH 25 Jul 2003 */
01159    char *sexi_end;
01160 
01161 ENTRY("mri_imcount_dicom") ;
01162 
01163    if( str_sexinfo != NULL ){ free(str_sexinfo); str_sexinfo=NULL; }
01164 
01165    if( !mri_possibly_dicom(fname) ) RETURN(0) ;  /* 07 May 2003 */
01166 
01167    /* extract the header from the file (cf. mri_dicom_hdr.[ch]) */
01168 
01169    mri_dicom_nohex(1) ; mri_dicom_noname(1) ;
01170    ppp = mri_dicom_header( fname ) ;
01171    if( ppp == NULL ) RETURN(0);
01172 
01173    /* find out where the pixel array is in the file */
01174 
01175    mri_dicom_pxlarr( &poff , &plen ) ;
01176    if( plen <= 0 ){ free(ppp) ; RETURN(0); }
01177 
01178    /* check if file is actually this big */
01179 
01180    { unsigned int psiz , fsiz ;
01181      fsiz = THD_filesize( fname ) ;
01182      psiz = (unsigned int)(poff) + plen ;
01183      if( fsiz < psiz ){ free(ppp) ; RETURN(0); }
01184    }
01185 
01186    /* find positions in header of elements we care about */
01187 
01188    for( ee=0 ; ee < NUM_ELIST ; ee++ )
01189      epos[ee] = strstr(ppp,elist[ee]) ;
01190 
01191    /* see if the header has the elements we absolutely need */
01192 
01193    if( epos[E_ROWS]           == NULL ||
01194        epos[E_COLUMNS]        == NULL ||
01195        epos[E_BITS_ALLOCATED] == NULL   ){ free(ppp) ; RETURN(0); }
01196 
01197    /* check if we have 1 sample per pixel (can't deal with 3 or 4 now) */
01198 
01199    if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
01200       ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
01201       if( ddd == NULL ){ free(ppp) ; RETURN(0); }
01202       ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
01203       if( ii != 1 ){ free(ppp) ; RETURN(0); }
01204    }
01205 
01206    /* check if photometric interpretation is MONOCHROME (don't like PALETTE) */
01207 
01208    if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
01209       ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
01210       if( ddd == NULL ){ free(ppp) ; RETURN(0); }
01211    }
01212 
01213    /* check if we have 8, 16, or 32 bits per pixel */
01214 
01215    ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
01216    if( ddd == NULL ){ free(ppp); RETURN(0); }
01217    bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
01218    switch( bpp ){
01219       default: free(ppp) ; RETURN(0);   /* bad value */
01220       case  8: datum = MRI_byte ; break ;
01221       case 16: datum = MRI_short; break ;
01222       case 32: datum = MRI_int  ; break ;
01223    }
01224    bpp /= 8 ; /* now bytes per pixel, instead of bits */
01225 
01226    /* get nx, ny, nz */
01227 
01228    ddd = strstr(epos[E_ROWS],"//") ;
01229    if( ddd == NULL ){ free(ppp) ; RETURN(0); }
01230    nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
01231    if( nx < 2 ){ free(ppp) ; RETURN(0); }
01232 
01233    ddd = strstr(epos[E_COLUMNS],"//") ;
01234    if( ddd == NULL ){ free(ppp) ; RETURN(0); }
01235    ny = 0 ; sscanf(ddd+2,"%d",&ny) ;
01236    if( ny < 2 ){ free(ppp) ; RETURN(0); }
01237 
01238    nz = 0 ;
01239    if( epos[E_NUMBER_OF_FRAMES] != NULL ){
01240      ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
01241      if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
01242    }
01243    if( nz == 0 ) nz = plen / (bpp*nx*ny) ;
01244 
01245    /*-- 28 Oct 2002: Check if this is a Siemens mosaic.        --*/
01246    /*-- 02 Dec 2002: Don't use Acquisition Matrix anymore;
01247                      instead, use the Siemens extra info
01248                      in epos[E_SIEMENS_2].                     --*/
01249    /*-- 24 Dec 2002: Extract str_sexinfo even if not a mosaic. --*/
01250 
01251    if(        epos[E_ID_MANUFACTURER]            != NULL &&
01252        strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
01253               epos[E_SIEMENS_2]                  != NULL    ){
01254 
01255      int len=0,loc=0 , aa,bb ;
01256      sscanf(epos[E_SIEMENS_2],"%x%x%d [%d" , &aa,&bb , &len,&loc ) ;
01257      if( len > 0 && loc > 0 ){
01258        FILE *fp = fopen( fname , "rb" ) ;
01259        if( fp != NULL ){
01260          str_sexinfo = extract_bytes_from_file( fp, (off_t)loc, (size_t)len, 1 ) ;
01261          fclose(fp) ;
01262        }
01263      }
01264    }
01265 
01266    if(        epos[E_ID_IMAGE_TYPE]              != NULL &&
01267        strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC")    != NULL &&
01268        str_sexinfo                               != NULL   ){
01269 
01270      /* 31 Oct 2002: extract extra Siemens info from file */
01271 
01272      /* KRH 25 Jul 2003 if start and end markers are present for
01273       * Siemens extra info, cut string down to those boundaries */
01274 
01275      sexi_start = strstr(str_sexinfo, "### ASCCONV BEGIN ###");
01276      sexi_end = strstr(str_sexinfo, "### ASCCONV END ###");
01277      if ((sexi_start != NULL) && (sexi_end != NULL)) {
01278        char *sexi_tmp;
01279        int sexi_size;
01280 
01281        sexi_size = sexi_end - sexi_start + 19 ;
01282        sexi_tmp = AFMALL( char,  sexi_size );
01283        memcpy(sexi_tmp,sexi_start,sexi_size);
01284        free(str_sexinfo);
01285        str_sexinfo = sexi_tmp;
01286      }
01287      /* end KRH 25 Jul 2003 change */
01288 
01289      sexinfo.good = 0 ;  /* start by marking it as bad */
01290      get_siemens_extra_info( str_sexinfo , &sexinfo ) ;
01291 
01292      if( sexinfo.good ){                                   /* if data is good */
01293 
01294        nz = sexinfo.mosaic_num ;
01295 
01296      } /* end of if sexinfo was good */
01297 
01298      else {                  /* warn if sexinfo was bad */
01299        static int nwarn=0 ;
01300        if( nwarn < NWMAX )
01301          fprintf(stderr,"++ DICOM WARNING: indecipherable SIEMENS MOSAIC info (%s) in file %s\n",
01302                  elist[E_SIEMENS_2] , fname ) ;
01303        if( nwarn == NWMAX )
01304          fprintf(stderr,"++ DICOM NOTICE: no more SIEMENS MOSAIC info messages will be printed\n");
01305        nwarn++ ;
01306      }
01307 
01308    } /* end of if a Siemens mosaic */
01309 
01310    free(ppp) ; RETURN(nz);
01311 }
01312 
01313 /*--------------------------------------------------------------------------------*/
01314 /*! Read some bytes from an open file at a given offset.  Return them in a
01315     newly malloc()-ed array.  If return value is NULL, something bad happened.
01316 ----------------------------------------------------------------------------------*/
01317 
01318 static char * extract_bytes_from_file( FILE *fp, off_t start, size_t len, int strize )
01319 {
01320    char *ar ;
01321    size_t nn , ii ;
01322 
01323    if( fp == NULL || len == 0 ) return NULL ;    /* bad inputs? */
01324    ar = AFMALL(char, len+1) ;                   /* make space for data */
01325    lseek( fileno(fp) , start , SEEK_SET ) ;      /* set file position */
01326    nn = fread( ar , 1 , len , fp ) ;             /* read data */
01327    if( nn == 0 ){ free(ar); return NULL; }       /* bad read? */
01328 
01329    if( strize ){                                 /* convert to C string? */
01330      for( ii=0 ; ii < nn ; ii++ )                /* scan for NULs and    */
01331        if( ar[ii] == '\0' ) ar[ii] = ' ' ;       /* replace with blanks  */
01332    }
01333    return ar ;
01334 }
01335 
01336 /*--------------------------------------------------------------------------------*/
01337 /*! Parse the Siemens extra stuff for mosaic information.
01338     Ad hoc, based on sample data and no documentation.
01339 ----------------------------------------------------------------------------------*/
01340 
01341 static void get_siemens_extra_info( char *str , Siemens_extra_info *mi )
01342 {
01343    char *cpt , *dpt ;
01344    int nn , mm , snum , last_snum=-1 ;
01345    int have_x[2] = {0,0},
01346        have_y[2] = {0,0},
01347        have_z[2] = {0,0};
01348    float x,y,z , val ;
01349    char name[1024] ;
01350 
01351    /*-- check for good inputs --*/
01352 
01353    if( mi == NULL ) return ;
01354 
01355    mi->good = 0 ;
01356    for( snum=0 ; snum < NMOMAX ; snum++ )
01357      mi->slice_xyz[snum][0] = mi->slice_xyz[snum][1] = mi->slice_xyz[snum][2] = 0.0 ;
01358 
01359    if( str == NULL || *str == '\0' ) return ;
01360 
01361    /*-- find string that starts the slice information array --*/
01362    /* 04 Mar 2003 reworked this section to skip "fake matches" *
01363     * of the target string present in some mosaic files in the *
01364     * binary section                                     --KRH */
01365    nn = 0;
01366    while (nn == 0) {
01367 
01368      cpt = strstr( str , "sSliceArray.asSlice[" ) ;
01369      if( cpt == NULL ) return ;
01370      /* interepret next string into
01371          snum = slice subscript (0,1,...)
01372          name = variable name
01373          val  = number of RHS of '=' sign
01374          mm   = # of bytes used in scanning the above */
01375 
01376      nn = sscanf( cpt , "sSliceArray.asSlice[%d].%1022s =%f%n" ,
01377                   &snum , name , &val , &mm ) ;
01378      str += 20; /* skip to end of "false match" string KRH */
01379    }
01380 
01381 
01382    /*-- scan for coordinates, until can't find a good string to scan --*/
01383 
01384 #if 0
01385    /* 25 Feb 2003 Changed logic here KRH */
01386    have_x = have_y = have_z = 0 ;
01387 #endif
01388 
01389    while(1){
01390 
01391 
01392 
01393      if( nn   <  3                   ) break ;  /* bad conversion set */
01394      if( snum <  0 || snum >= NMOMAX ) break ;  /* slice number out of range */
01395 
01396 /* 21 Feb 2003 rework this section to allow for missing coordinates in
01397  * some mosaic files                                        --KRH  */
01398 #if 0
01399      /* assign val based on name */
01400 
01401           if( strcmp(name,"sPosition.dSag") == 0 ){ x = val; have_x = 1; }
01402      else if( strcmp(name,"sPosition.dCor") == 0 ){ y = val; have_y = 1; }
01403      else if( strcmp(name,"sPosition.dTra") == 0 ){ z = val; have_z = 1; }
01404 
01405      /* if now have all 3 coordinates, save them */
01406 
01407      if( have_x && have_y && have_z ){
01408        mi->slice_xyz[snum][0] = x; mi->slice_xyz[snum][1] = y; mi->slice_xyz[snum][2] = z;
01409        last_snum = snum ;
01410        have_x = have_y = have_z = 0 ;
01411      }
01412 #endif
01413 
01414      if( strcmp(name,"sPosition.dSag") == 0 ){ 
01415        mi->slice_xyz[snum][0] = val;
01416        if (snum < 2) have_x[snum] = 1;
01417        last_snum = snum;
01418      } else if( strcmp(name,"sPosition.dCor") == 0 ){ 
01419        mi->slice_xyz[snum][1] = val;
01420        if (snum < 2) have_y[snum] = 1;
01421        last_snum = snum;
01422      } else if( strcmp(name,"sPosition.dTra") == 0 ){ 
01423        mi->slice_xyz[snum][2] = val;
01424        if (snum < 2) have_z[snum] = 1;
01425        last_snum = snum;
01426      }
01427 
01428      /* skip to next slice array assignment string (which may not be a coordinate) */
01429 
01430      dpt = cpt + mm ;                                           /* just after 'val' */
01431      cpt =  dpt ;
01432      while( isspace(*cpt) ) cpt++ ;                             /* skip over whitespace */
01433      if( cpt-dpt > 16 ) break ;                                 /* too much space */
01434      if( strncmp(cpt,"sSliceArray.asSlice[",20) != 0 ) break ;   /* bad next line */
01435      /* 04 Mar 2003 moved this stuff around to allow for locating "fake matches"  *
01436       * of the target text in some mosaic files' binary sections                  */
01437 
01438      /* interpret next string into
01439          snum = slice subscript (0,1,...)
01440          name = variable name
01441          val  = number of RHS of '=' sign
01442          mm   = # of bytes used in scanning the above */
01443 
01444      nn = sscanf( cpt , "sSliceArray.asSlice[%d].%1022s =%f%n" ,
01445                   &snum , name , &val , &mm ) ;
01446 
01447    }
01448 
01449    /* if got at least 1 slice info, mark data as being good */
01450 
01451    if( last_snum >= 0 ){
01452      mi->good       = 1 ;
01453      if (have_x[0] && have_x[1]) mi->have_data[0] = 1;
01454      if (have_y[0] && have_y[1]) mi->have_data[1] = 1;
01455      if (have_z[0] && have_z[1]) mi->have_data[2] = 1;
01456      mi->mosaic_num = last_snum+1 ;
01457    }
01458 
01459    return ;
01460 }
01461 
01462 /*--------------------------------------------------------------------------*/
01463 /*! Test if a file is possibly a DICOM file.  -- RWCox - 07 May 2003
01464 ----------------------------------------------------------------------------*/
01465 
01466 int mri_possibly_dicom( char *fname )
01467 {
01468 #undef  BSIZ
01469 #define BSIZ 4096
01470    FILE *fp ;
01471    unsigned char buf[BSIZ] , *cpt ;
01472    int nn , ii ;
01473 
01474    if( fname == NULL || *fname == '\0' ) return 0 ;
01475    fp = fopen( fname , "rb" ) ; if( fp == NULL ) return 0 ;
01476 
01477    /* read 1st buffer */
01478 
01479    nn = fread( buf , 1 , BSIZ , fp ) ;
01480    if( nn < 256 ){ fclose(fp) ; return 0 ; }  /* too short */
01481 
01482    /* easy: check if has 'DICM' marker at offset 128..131 */
01483 
01484    if( buf[128]=='D' && buf[129]=='I' && buf[130]=='C' && buf[131]=='M' ){
01485      fclose(fp) ; return 1 ;
01486    }
01487 
01488    /* hard: scan file for sequence: E0 7F 10 00 (image data attribute) */
01489 
01490    while(1){
01491 
01492      cpt = memchr( buf, 0xe0, nn ) ;                /* look for E0 */
01493 
01494      if( cpt == NULL ){                        /* skip this buffer */
01495        nn = fread( buf , 1 , BSIZ , fp ) ;      /* and get another */
01496        if( nn < 256 ){ fclose(fp) ; return 0 ; }
01497        continue ;
01498      }
01499 
01500      ii = nn - (cpt-buf) ;               /* num char to end of buf */
01501      if( ii <= 4 ){                     /* too close to end of buf */
01502        memmove( buf , cpt , ii ) ;
01503        nn = fread( buf+ii , 1 , BSIZ-ii , fp ) ; nn += ii ;
01504        if( nn < 256 ){ fclose(fp) ; return 0 ; }
01505        cpt = buf ; ii = nn ;
01506      }
01507 
01508      /* see if we got what we want */
01509 
01510      if( *cpt==0xe0 && *(cpt+1)==0x7f && *(cpt+2)==0x10 && *(cpt+3)==0x00 ){
01511        fclose(fp) ; return 1 ;
01512      }
01513 
01514      /* no?  start again at next char in buf */
01515 
01516      memmove( buf , cpt+1 , ii-1 ) ; nn = ii-1 ;
01517    }
01518 }