02743 {
02744     static struct stat buf ; int ii ;
02745 
02746     if( pathname == NULL ) return -1 ;
02747         ii = stat( pathname , &buf ) ; if( ii != 0 ) return -1 ;
02748 
02749     return buf.st_size ;
02750 }

00205 {
00206     DCM_OBJECT * object;
00207     CONDITION cond;
00208     CTNBOOLEAN verbose = FALSE ,
00209                exitFlag = FALSE,
00210                formatFlag = FALSE;
00211     unsigned long
00212         options = DCM_ORDERLITTLEENDIAN;
00213     long vmLimit = rwc_vm ;             /* 28 Oct 2002 */
00214     LST_HEAD* fileNames = 0;
00215     UTL_FILEITEM* p = NULL;
00216 
00217     char *ppp=NULL ;
00218 
00219 ENTRY("mri_dicom_header") ;
00220 
00221     if( fname == NULL ) RETURN(NULL) ;
00222 
00223     RWC_set_endianosity() ;
00224 
00225     { char *eee = getenv("AFNI_TRACE") ;
00226       if( eee!=NULL && (*eee=='y' || *eee=='Y') ) verbose = TRUE ;
00227     }
00228 
00229     DCM_Debug(verbose);
00230 
00231     RWC_clear_pbuf() ; pxl_len = 0 ; pxl_off = 0 ;
00232 
00233 STATUS(fname) ;
00234     rwc_fd = -1 ;
00235     cond = DCM_OpenFile(fname, options, &object);
00236     if (cond != DCM_NORMAL && ((options & DCM_PART10FILE) == 0)) {
00237 STATUS("DCM_OpenFile open failed; try again as Part 10") ;
00238       (void) DCM_CloseObject(&object);
00239       (void) COND_PopCondition(TRUE);
00240       if( rwc_fd >= 0 ){ close(rwc_fd); rwc_fd = -1; }
00241       cond = DCM_OpenFile(fname, options | DCM_PART10FILE, &object);
00242     }
00243     if (cond == DCM_NORMAL) {
00244 STATUS("DCM_OpenFile is good") ;
00245        RWC_printf("DICOM File: %s\n", fname);
00246        if (formatFlag)
00247          cond = DCM_FormatElements(&object, vmLimit, "");
00248        else
00249          cond = DCM_DumpElements(&object, vmLimit);
00250     } else {
00251 STATUS("DCM_OpenFile failed") ;
00252     }
00253     (void) DCM_CloseObject(&object);
00254     (void) COND_PopCondition(TRUE);
00255 
00256     if( pbuf != NULL ){
00257       ppp = strdup(pbuf) ; RWC_clear_pbuf() ;
00258     }
00259 
00260     if( rwc_fd >= 0 ){ close(rwc_fd); rwc_fd = -1; }
00261 
00262     RETURN(ppp);
00263 }

00175 {
00176    if( ii )
00177      rwc_opt |= RWC_NOHEX_MASK ;
00178    else if( rwc_opt & RWC_NOHEX_MASK )
00179      rwc_opt ^= RWC_NOHEX_MASK ;
00180 }

00167 {
00168    if( ii )
00169      rwc_opt |= RWC_NONAME_MASK ;
00170    else if( rwc_opt & RWC_NONAME_MASK )
00171      rwc_opt ^= RWC_NONAME_MASK ;
00172 }

00155 {
00156    *poff = pxl_off ; *plen = pxl_len ;
00157 }

00208 {
00209 #undef  BSIZ
00210 #define BSIZ 4096
00211    FILE *fp ;
00212    unsigned char buf[BSIZ] , *cpt ;
00213    int nn , ii ;
00214 
00215    if( fname == NULL || *fname == '\0' ) return 0 ;
00216    fp = fopen( fname , "rb" ) ; if( fp == NULL ) return 0 ;
00217 
00218    /* read 1st buffer */
00219 
00220    nn = fread( buf , 1 , BSIZ , fp ) ;
00221    if( nn < 256 ){ fclose(fp) ; return 0 ; }  /* too short */
00222 
00223    /* easy: check if has 'DICM' marker at offset 128..131 */
00224 
00225    if( buf[128]=='D' && buf[129]=='I' && buf[130]=='C' && buf[131]=='M' ){
00226      fclose(fp) ; return 1 ;
00227    }
00228 
00229    /* hard: scan file for sequence: E0 7F 10 00 (image data attribute) */
00230 
00231    while(1){
00232 
00233      cpt = memchr( buf, 0xe0, nn ) ;                /* look for E0 */
00234 
00235      if( cpt == NULL ){                        /* skip this buffer */
00236        nn = fread( buf , 1 , BSIZ , fp ) ;      /* and get another */
00237        if( nn < 256 ){ fclose(fp) ; return 0 ; }
00238        continue ;
00239      }
00240 
00241      ii = nn - (cpt-buf) ;               /* num char to end of buf */
00242      if( ii <= 4 ){                     /* too close to end of buf */
00243        memmove( buf , cpt , ii ) ;
00244        nn = fread( buf+ii , 1 , BSIZ-ii , fp ) ; nn += ii ;
00245        if( nn < 256 ){ fclose(fp) ; return 0 ; }
00246        cpt = buf ; ii = nn ;
00247      }
00248 
00249      /* see if we got what we want */
00250 
00251      if( *cpt==0xe0 && *(cpt+1)==0x7f && *(cpt+2)==0x10 && *(cpt+3)==0x00 ){
00252        fclose(fp) ; return 1 ;
00253      }
00254 
00255      /* no?  start again at next char in buf */
00256 
00257      memmove( buf , cpt+1 , ii-1 ) ; nn = ii-1 ;
00258    }
00259 }

00265 {
00266    char *ppp , *ddd ;
00267    off_t poff ;
00268    unsigned int plen ;
00269    char *epos[NUM_ELIST] ;
00270    int ii,jj , ee , bpp , datum ;
00271    int nx,ny,nz , swap ;
00272    float dx,dy,dz,dt ;
00273    MRI_IMAGE *im ;
00274    void *iar ;
00275    FILE *fp ;
00276    int have_orients=0 ;
00277    int ior,jor,kor ;
00278    static int nzoff=0 ;   /* for determining z-axis orientation/offset from multiple calls */
00279 
00280    float dxx,dyy,dzz ;
00281 
00282    char *eee ;
00283    float rescale_slope=0.0 , rescale_inter=0.0 ;  /* 23 Dec 2002 */
00284 
00285    ENTRY("mri_read_dicom") ;
00286 
00287    if( !mri_possibly_dicom(fname) )
00288    {
00289       if(debug > 2) fprintf(stderr,"** MRD: mri_possibly_dicom() failure\n");
00290       RETURN(NULL) ;  /* 07 May 2003 */
00291    }
00292 
00293    /* extract header info from file into a string
00294       - cf. mri_dicom_hdr.[ch]
00295       - run 'dicom_hdr -noname fname' to see the string format */
00296 
00297    mri_dicom_nohex(1) ;              /* don't print ints in hex mode */
00298    mri_dicom_noname(1) ;             /* don't print names, just tags */
00299    ppp = mri_dicom_header( fname ) ; /* print header to malloc()-ed string */
00300    if( ppp == NULL )
00301    {
00302       if(debug > 2) fprintf(stderr,"** MRD: mri_dicom_hdr() failure\n");
00303       RETURN(NULL);   /* didn't work; not a DICOM file? */
00304    }
00305 
00306    /* find out where the pixel array is in the file */
00307 
00308    mri_dicom_pxlarr( &poff , &plen ) ;
00309    if( plen <= 0 ){
00310       if(debug > 2) fprintf(stderr,"** MRD: bad plen, %u\n", plen);
00311       free(ppp) ;
00312       RETURN(NULL);
00313    }
00314    if( debug > 3 )
00315       fprintf(stderr,"-d dicom: poff, plen = %d, %d\n",(int)poff,(int)plen);
00316 
00317    /* check if file is actually this big (maybe it was truncated) */
00318 
00319    { unsigned int psiz , fsiz ;
00320      fsiz = l_THD_filesize( fname ) ;
00321      psiz = (unsigned int)(poff) + plen ;
00322      if( fsiz < psiz ){
00323         if(debug > 2) fprintf(stderr,"** MRD: fsiz < psiz (%d,%d)\n",fsiz,psiz);
00324         free(ppp) ;
00325         RETURN(NULL);
00326      }
00327    }
00328 
00329    /* find positions in header of elements we care about */
00330 
00331    for( ee=0 ; ee < NUM_ELIST ; ee++ )
00332      epos[ee] = strstr(ppp,elist[ee]) ;
00333 
00334    /* see if the header has the elements we absolutely need */
00335 
00336    if( epos[E_ROWS]           == NULL ||
00337        epos[E_COLUMNS]        == NULL ||
00338        epos[E_BITS_ALLOCATED] == NULL   ){
00339       if(debug > 2)
00340          fprintf(stderr,"** MRD: missing ROWS, COLS or BITS (%p,%p,%p)\n",
00341                  epos[E_ROWS], epos[E_COLUMNS], epos[E_BITS_ALLOCATED]);
00342       free(ppp) ;
00343       RETURN(NULL);
00344    }
00345 
00346    /* check if we have 1 sample per pixel (can't deal with 3 or 4 now) */
00347 
00348    if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
00349       ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
00350       if( ddd == NULL ){
00351          if(debug > 2) fprintf(stderr,"** MRD: missing E_SAMPLES_PER_PIXEL\n");
00352          free(ppp) ;
00353          RETURN(NULL);
00354       }
00355       ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
00356        if( ii != 1 ){
00357          if(debug > 2) fprintf(stderr,"** MRD: SAM_PER_PIX != 1, %d\n",ii);
00358          free(ppp) ;
00359          RETURN(NULL);
00360       }
00361    }
00362 
00363    /* check if photometric interpretation is MONOCHROME (don't like PALETTE) */
00364 
00365    if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
00366       ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
00367       if( ddd == NULL ){
00368       if(debug > 2) fprintf(stderr,"** MRD: photometric not monochrome\n");
00369         free(ppp) ;
00370         RETURN(NULL);
00371       }
00372    }
00373 
00374    /* check if we have 8, 16, or 32 bits per pixel */
00375 
00376    ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
00377    if( ddd == NULL ){
00378       if(debug > 2) fprintf(stderr,"** MRD: missing BITS_ALLOCATED\n");
00379       free(ppp);
00380       RETURN(NULL);
00381    }
00382    bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
00383    switch( bpp ){
00384       default:
00385         if(debug > 2) fprintf(stderr,"** MRD: bad bpp value, %d\n",bpp);
00386         free(ppp); RETURN(NULL);    /* bad value */
00387       case  8: datum = MRI_byte ; break ;
00388       case 16: datum = MRI_short; break ;
00389       case 32: datum = MRI_int  ; break ;  /* probably not present in DICOM? */
00390    }
00391    bpp /= 8 ; /* now bytes per pixel, instead of bits */
00392 
00393    if( debug > 3 ) fprintf(stderr,"-d dicom: datum %d\n",datum);
00394 
00395    /*** Print some warnings if appropriate ***/
00396 
00397    /* check if BITS_STORED and HIGH_BIT are aligned */
00398 
00399 #define NWMAX 2
00400    if( epos[E_BITS_STORED] != NULL && epos[E_HIGH_BIT] != NULL ){
00401      int bs=0 , hb=0 ;
00402      ddd = strstr(epos[E_BITS_STORED],"//") ; sscanf(ddd+2,"%d",&bs) ;
00403      ddd = strstr(epos[E_HIGH_BIT],"//")    ; sscanf(ddd+2,"%d",&hb) ;
00404      if( bs != hb+1 ){
00405        static int nwarn=0 ;
00406        if( nwarn < NWMAX )
00407          fprintf(stderr,
00408                  "++ DICOM WARNING: file %s has Bits_Stored=%d and High_Bit=%d\n",
00409                  fname,bs,hb) ;
00410        if( nwarn == NWMAX )
00411          fprintf(stderr,"++ DICOM WARNING: no more Bits_Stored messages will be printed\n") ;
00412        nwarn++ ;
00413      }
00414    }
00415 
00416    /* check if Rescale is ordered */
00417    /* 23 Dec 2002: actually get the rescale params, if environment says to */
00418 
00419    if( epos[E_RESCALE_INTERCEPT] != NULL && epos[E_RESCALE_SLOPE] != NULL ){
00420      if( debug > 0 ){
00421        static int nwarn=0 ;
00422        if( nwarn == 0 ){
00423          fprintf(stderr,"++ DICOM file has Rescale tags, enforcing them...\n");
00424          fprintf(stderr,"   (no more Rescale tags messages will be printed)\n");
00425        }
00426        nwarn++ ;
00427      }
00428      ddd = strstr(epos[E_RESCALE_INTERCEPT],"//") ;
00429      sscanf(ddd+2,"%f",&rescale_inter) ;
00430      ddd = strstr(epos[E_RESCALE_SLOPE    ],"//") ;
00431      sscanf(ddd+2,"%f",&rescale_slope) ;
00432    }
00433 
00434    /*** extract attributes of the image(s) to be read in ***/
00435 
00436    /* get image nx & ny */
00437 
00438    ddd = strstr(epos[E_ROWS],"//") ;                         /* 31 Oct 2002: */
00439    if( ddd == NULL ){                 /* Oops: ROWS is ny and COLUMNS is nx! */
00440       if(debug > 2) fprintf(stderr,"** MRD: missing E_ROWS\n");
00441       free(ppp) ;
00442       RETURN(NULL);
00443    }
00444    ny = 0 ; sscanf(ddd+2,"%d",&ny) ;
00445    if( ny < 2 ){
00446       if(debug > 2) fprintf(stderr,"** MRD: bad ny = %d\n",ny);
00447       free(ppp) ;
00448       RETURN(NULL);
00449    }
00450 
00451    ddd = strstr(epos[E_COLUMNS],"//") ;
00452    if( ddd == NULL ){
00453       if(debug > 2) fprintf(stderr,"** MRD: missing E_COLUMNS\n");
00454       free(ppp) ;
00455       RETURN(NULL);
00456    }
00457    nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
00458    if( nx < 2 ){
00459       if(debug > 2) fprintf(stderr,"** MRD: bad nx = %d\n",nx);
00460       free(ppp) ;
00461       RETURN(NULL);
00462    }
00463 
00464    /* get number of slices */
00465 
00466    nz = 0 ;
00467    if( epos[E_NUMBER_OF_FRAMES] != NULL ){
00468      ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
00469      if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
00470      if(debug>3) fprintf(stderr,"-d number of frames = %d\n",nz);
00471    }
00472 
00473    /* if didn't get nz above, make up a value */
00474 
00475    if( nz == 0 ) nz = plen / (bpp*nx*ny) ;   /* compute from image array size */
00476    if( nz == 0 ){
00477       if(debug > 2) fprintf(stderr,"** MRD: bad nz = %d\n", nz);
00478       free(ppp) ;
00479       RETURN(NULL);
00480    }
00481 
00482    if( nz != 1 ){
00483       fprintf(stderr,"** MRD: nz = %d, plen,bpp,nx,ny = %d,%d,%d,%d\n",
00484          nz, plen, bpp, nx, ny);
00485       free(ppp);
00486       RETURN(NULL);
00487    }
00488 
00489    /*-- 28 Oct 2002: Check if this is a Siemens mosaic.        --*/
00490    /*-- 02 Dec 2002: Don't use Acquisition Matrix anymore;
00491                      instead, use the Siemens extra info
00492                      in epos[E_SIEMENS_2].                     --*/
00493    /*-- 24 Dec 2002: Extract str_sexinfo even if not a mosaic. --*/
00494 
00495    /* ignore sexinfo for Dimon, since it becomes a mosaic-only case */
00496 
00497    if(        epos[E_ID_MANUFACTURER]            != NULL &&
00498        strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
00499               epos[E_SIEMENS_2]                  != NULL    ){
00500       if( debug > 3 ) fprintf(stderr,"-d ignoring sexinfo\n");
00501    }
00502 
00503    /*-- try to get dx, dy, dz, dt --*/
00504 
00505    dx = dy = dz = dt = 0.0 ;
00506 
00507    /* dx,dy first */
00508 
00509    if( epos[E_PIXEL_SPACING] != NULL ){
00510      ddd = strstr(epos[E_PIXEL_SPACING],"//") ;
00511      if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00512      if( debug > 3 )
00513         fprintf(stderr,"-d dicom: read dx,dy from PIXEL_SP: %f, %f\n", dx, dy);
00514      if( dy == 0.0 && dx > 0.0 ) dy = dx ;
00515    }
00516    if( dx == 0.0 && epos[E_FIELD_OF_VIEW] != NULL ){
00517      ddd = strstr(epos[E_FIELD_OF_VIEW],"//") ;
00518      if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00519      if( debug > 3 )
00520         fprintf(stderr,"-d dicom: read dx,dy from FOV: %f, %f\n", dx, dy);
00521      if( dx > 0.0 ){
00522        if( dy == 0.0 ) dy = dx ;
00523        dx /= nx ; dy /= ny ;
00524      }
00525    }
00526 
00527    /*-- 27 Nov 2002: fix stupid GE error,
00528                      where the slice spacing is really the slice gap --*/
00529 
00530    { int stupid_ge_fix , no_stupidity ;
00531      float sp=0.0 , th=0.0 ;
00532      static int nwarn=0 ;
00533                                                            /* 03 Mar 2003 */
00534      eee           = getenv("AFNI_SLICE_SPACING_IS_GAP") ;
00535      stupid_ge_fix = (eee != NULL && (*eee=='Y' || *eee=='y') );
00536      no_stupidity  = (eee != NULL && (*eee=='N' || *eee=='n') );
00537 
00538      if( epos[E_SLICE_SPACING] != NULL ){    /* get reported slice spacing */
00539        ddd = strstr(epos[E_SLICE_SPACING],"//") ;
00540        if( ddd != NULL ) sscanf( ddd+2 , "%f" , &sp ) ;
00541      }
00542      if( epos[E_SLICE_THICKNESS] != NULL ){  /* get reported slice thickness */
00543        ddd = strstr(epos[E_SLICE_THICKNESS],"//") ;
00544        if( ddd != NULL ) sscanf( ddd+2 , "%f" , &th ) ;
00545      }
00546 
00547      if(debug>3) fprintf(stderr,"-d dicom: thick, spacing = %f, %f\n", th, sp);
00548 
00549      th = fabs(th) ; sp = fabs(sp) ;         /* we don't use the sign */
00550 
00551      if( stupid_ge_fix ){                    /* always be stupid */
00552        dz = sp+th ;
00553      } else {
00554 
00555        if( no_stupidity && sp > 0.0 )        /* 13 Jan 2004: if 'NO', then */
00556          dz = sp ;                           /* always use spacing if present */
00557        else
00558          dz = (sp > th) ? sp : th ;          /* the correct-ish DICOM way */
00559 
00560 #define GFAC 0.99
00561 
00562        if( !no_stupidity ){                 /* unless stupidity is turned off */
00563          if( sp > 0.0 && sp < GFAC*th ) dz = sp+th;/* the stupid GE way again */
00564 
00565          if( sp > 0.0 && sp < GFAC*th && nwarn < NWMAX ){
00566            fprintf(stderr,
00567                    "++ DICOM WARNING: file %s has Slice_Spacing=%f smaller than Slice_Thickness=%f\n",
00568                    fname , sp , th ) ;
00569            if( nwarn == 0 )
00570             fprintf(stderr,
00571               "\n"
00572       "++  Setting environment variable AFNI_SLICE_SPACING_IS_GAP       ++\n"
00573       "++   to YES will make the center-to-center slice distance        ++\n"
00574       "++   be set to Slice_Spacing+Slice_Thickness=%6.3f.             ++\n"
00575       "++  This is against the DICOM standard [attribute (0018,0088)    ++\n"
00576       "++   is defined as the center-to-center spacing between slices,  ++\n"
00577       "++   NOT as the edge-to-edge gap between slices], but it seems   ++\n"
00578       "++   to be necessary for some GE scanners.                       ++\n"
00579       "++                                                               ++\n"
00580       "++  This correction has been made on this data: dz=%6.3f.       ++\n"
00581       "++                                                               ++\n"
00582       "++  Setting AFNI_SLICE_SPACING_IS_GAP to NO means that the       ++\n"
00583       "++  DICOM Slice_Spacing variable will be used for dz, replacing  ++\n"
00584       "++  the Slice_Thickness variable.  This usage may be required    ++\n"
00585       "++  for some pulse sequences on Phillips scanners.               ++\n"
00586               "\n\a" ,
00587              sp+th , dz ) ;
00588          }
00589          if( sp > 0.0 && sp < th && nwarn == NWMAX )
00590            fprintf(stderr,"++ DICOM WARNING: no more Slice_Spacing messages will be printed\n") ;
00591          nwarn++ ;
00592        }
00593      }
00594      if( dz == 0.0 && dx != 0.0 ) dz = 1.0 ;               /* nominal dz */
00595 
00596    } /*-- end of dz code, with all its stupidities --*/
00597 
00598    if(debug > 3) fprintf(stderr,"-d dicom: using dxyz = %f, %f, %f\n",dx,dy,dz);
00599 
00600    /* get dt */
00601 
00602    if( epos[E_REPETITION_TIME] != NULL ){
00603      ddd = strstr(epos[E_REPETITION_TIME],"//") ;
00604      if( ddd != NULL ){
00605        sscanf( ddd+2 , "%f" , &dt ) ;
00606        dt *= 0.001 ;   /* ms to s */
00607        if(debug > 3) fprintf(stderr,"-d dicom: rep. time dt = %f sec.\n",dt);
00608      }
00609    }
00610 
00611    /* check if we might have 16 bit unsigned data that fills all bits */
00612 
00613    if( epos[E_RS_STUDY_NUM] != NULL ){
00614      ddd = strstr(epos[E_RS_STUDY_NUM],"//") ;
00615      if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.study);
00616    }
00617 
00618    if( epos[E_RS_SERIES_NUM] != NULL ){
00619      ddd = strstr(epos[E_RS_SERIES_NUM],"//") ;
00620      if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.series);
00621    }
00622 
00623    if( epos[E_RS_IMAGE_NUM] != NULL ){
00624      ddd = strstr(epos[E_RS_IMAGE_NUM],"//") ;
00625      if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.image);
00626    }
00627 
00628    /** Finally! Read images from file. **/
00629 
00630    fp = fopen( fname , "rb" ) ;
00631    if( fp == NULL ){
00632       if(debug > 2) fprintf(stderr,"** MRD: failed to open file '%s'\n",fname);
00633       free(ppp) ;
00634       RETURN(NULL);
00635    }
00636    lseek( fileno(fp) , poff , SEEK_SET ) ;
00637 
00638    /* DICOM files are stored in LSB first (little endian) mode */
00639 
00640    set_endianosity() ; swap = !LITTLE_ENDIAN_ARCHITECTURE ;
00641 
00642    /* use pre-(28 Oct 2002) method, mosaics are not allowed now */
00643 
00644    { /* no mri_new(), as we don't want to link everything */
00645      im = (MRI_IMAGE *)calloc(1, sizeof(MRI_IMAGE));
00646      if( !im ){
00647         fprintf(stderr,"** MRD: im malloc failure!\n");
00648         free(ppp);
00649         RETURN(NULL);
00650      }
00651      im->nx = nx;  im->ny = ny;
00652      im->nxy = nx * ny;
00653      im->nz = im->nt = im->nu = im->nv = im->nw = 1;
00654      im->nxyz = im->nxyzt = im->nvox = im->nxy; 
00655      im->kind = datum;
00656      im->dx = im->dy = im->dz = im->dt = im->du = im->dv = 1.0;
00657      im->dw = -666.0;
00658      im->pixel_size = bpp;
00659      if( data ){ /* if data is not set, do not read it in */
00660         if( ! *data ){
00661            *data = calloc(im->nvox, im->pixel_size);
00662            if( debug > 3 )
00663               fprintf(stderr,"+d dicom: image data alloc: %d x %d bytes\n",
00664                  im->nvox, im->pixel_size);
00665         }
00666         if( ! *data ){
00667            fprintf(stderr,"** MRD: image data alloc failure\n");
00668            free(ppp);
00669            RETURN(NULL);
00670         } 
00671      }
00672    }
00673 
00674    if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
00675      im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
00676    }
00677    if( dt > 0.0 ) im->dt = dt ;
00678    
00679    if( !data ) fclose(fp) ; 
00680    else{
00681       iar = *data;
00682       fread( iar , bpp , nx*ny , fp ) ;    /* read data directly into it */
00683    
00684       if( swap ){                          /* swap bytes? */
00685         switch( im->pixel_size ){
00686           default: break ;
00687           case 2:  swap_twobytes (   im->nvox, iar ) ; break ;  /* short */
00688           case 4:  swap_fourbytes(   im->nvox, iar ) ; break ;  /* int, float */
00689           case 8:  swap_fourbytes( 2*im->nvox, iar ) ; break ;  /* complex */
00690         }
00691         im->was_swapped = 1 ;
00692       }
00693    
00694       /* store auxiliary data in image struct */
00695    
00696       fclose(fp) ;     /* 10 Sep 2002: oopsie - forgot to close file */
00697    
00698       /*-- 23 Dec 2002: implement Rescale, if ordered --*/
00699    
00700       if( rescale_slope > 0.0 ){
00701         for( ii=0 ; ii < 1 ; ii++ ){
00702           switch( im->kind ){
00703             case MRI_byte:{
00704               byte *ar = (byte *)*data;
00705               for( jj=0 ; jj < im->nvox ; jj++ )
00706                 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00707             }
00708             break ;
00709    
00710             case MRI_short:{
00711               short *ar = (short *)*data;
00712               for( jj=0 ; jj < im->nvox ; jj++ )
00713                 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00714             }
00715             break ;
00716    
00717             case MRI_int:{
00718               int *ar = (int *)*data;
00719               for( jj=0 ; jj < im->nvox ; jj++ )
00720                 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00721             }
00722             break ;
00723             default:{
00724               fprintf(stderr,"** MRD: bad kind case (%d) for rescale_slope\n",
00725                       im->kind);
00726               free(ppp);  free(*data);  *data = NULL;
00727               RETURN(NULL);
00728             }
00729             break ;
00730           }
00731         }
00732       } /* end of Rescale */
00733    }
00734 
00735    /*-- no WINDOW in Dimon */
00736 
00737    /*-- store some extra information in DI_MRL globals, too? --*/
00738 
00739    if( dt > 0.0 && DI_MRL_tr <= 0.0 ) DI_MRL_tr = dt ;
00740 
00741    /*-- try to get image orientation fields (also, set ior,jor,kor) --*/
00742 
00743    if( epos[E_IMAGE_ORIENTATION] != NULL ){
00744      /* direction cosines of image plane */
00745 
00746      ddd = strstr(epos[E_IMAGE_ORIENTATION],"//") ;
00747      if( ddd != NULL ){
00748        float xc1=0.0,xc2=0.0,xc3=0.0 , yc1=0.0,yc2=0.0,yc3=0.0 ;
00749        float xn,yn ; int qq ;
00750        qq=sscanf(ddd+2,"%f\\%f\\%f\\%f\\%f\\%f",&xc1,&xc2,&xc3,&yc1,&yc2,&yc3);
00751        xn = sqrt( xc1*xc1 + xc2*xc2 + xc3*xc3 ) ; /* vector norms */
00752        yn = sqrt( yc1*yc1 + yc2*yc2 + yc3*yc3 ) ;
00753        if( qq == 6 && xn > 0.0 && yn > 0.0 ){     /* both vectors OK */
00754 
00755          xc1 /= xn ; xc2 /= xn ; xc3 /= xn ;      /* normalize vectors */
00756          yc1 /= yn ; yc2 /= yn ; yc3 /= yn ;
00757 
00758          if( !use_DI_MRL_xcos ){
00759            DI_MRL_xcos[0] = xc1 ; DI_MRL_xcos[1] = xc2 ;  /* save direction */
00760            DI_MRL_xcos[2] = xc3 ; use_DI_MRL_xcos = 1 ;   /* cosine vectors */
00761          }
00762 
00763          if( !use_DI_MRL_ycos ){
00764            DI_MRL_ycos[0] = yc1 ; DI_MRL_ycos[1] = yc2 ;
00765            DI_MRL_ycos[2] = yc3 ; use_DI_MRL_ycos = 1 ;
00766          }
00767 
00768          /* x-axis orientation */
00769          /* ior determines which spatial direction is x-axis  */
00770          /* and is the direction that has the biggest change */
00771 
00772          dxx = fabs(xc1) ; ior = 1 ;
00773          dyy = fabs(xc2) ; if( dyy > dxx ){ ior=2; dxx=dyy; }
00774          dzz = fabs(xc3) ; if( dzz > dxx ){ ior=3;        }
00775          dxx = DI_MRL_xcos[ior-1] ; if( dxx < 0. ) ior = -ior;
00776 
00777          if( DI_MRL_orients[0] == '\0' ){
00778            switch( ior ){
00779              case -1: DI_MRL_orients[0] = 'L'; DI_MRL_orients[1] = 'R'; break;
00780              case  1: DI_MRL_orients[0] = 'R'; DI_MRL_orients[1] = 'L'; break;
00781              case -2: DI_MRL_orients[0] = 'P'; DI_MRL_orients[1] = 'A'; break;
00782              case  2: DI_MRL_orients[0] = 'A'; DI_MRL_orients[1] = 'P'; break;
00783              case  3: DI_MRL_orients[0] = 'I'; DI_MRL_orients[1] = 'S'; break;
00784              case -3: DI_MRL_orients[0] = 'S'; DI_MRL_orients[1] = 'I'; break;
00785              default: DI_MRL_orients[0] ='\0'; DI_MRL_orients[1] ='\0'; break;
00786            }
00787          }
00788 
00789          /* y-axis orientation */
00790          /* jor determines which spatial direction is y-axis  */
00791          /* and is the direction that has the biggest change */
00792 
00793          dxx = fabs(yc1) ; jor = 1 ;
00794          dyy = fabs(yc2) ; if( dyy > dxx ){ jor=2; dxx=dyy; }
00795          dzz = fabs(yc3) ; if( dzz > dxx ){ jor=3;        }
00796          dyy = DI_MRL_ycos[jor-1] ; if( dyy < 0. ) jor = -jor;
00797          if( DI_MRL_orients[2] == '\0' ){
00798            switch( jor ){
00799              case -1: DI_MRL_orients[2] = 'L'; DI_MRL_orients[3] = 'R'; break;
00800              case  1: DI_MRL_orients[2] = 'R'; DI_MRL_orients[3] = 'L'; break;
00801              case -2: DI_MRL_orients[2] = 'P'; DI_MRL_orients[3] = 'A'; break;
00802              case  2: DI_MRL_orients[2] = 'A'; DI_MRL_orients[3] = 'P'; break;
00803              case  3: DI_MRL_orients[2] = 'I'; DI_MRL_orients[3] = 'S'; break;
00804              case -3: DI_MRL_orients[2] = 'S'; DI_MRL_orients[3] = 'I'; break;
00805              default: DI_MRL_orients[2] ='\0'; DI_MRL_orients[3] ='\0'; break;
00806            }
00807          }
00808 
00809          DI_MRL_orients[6] = '\0' ;   /* terminate orientation string */
00810 
00811          kor = 6 - abs(ior)-abs(jor) ;   /* which spatial direction is z-axis */
00812                                          /* where 1=L-R, 2=P-A, 3=I-S */
00813          have_orients = 1 ;
00814 
00815          if( debug > 3 )
00816            fprintf(stderr,"-d dicom: DI_MRL_orients = %s, [ijk]or = %d,%d,%d\n",
00817                    DI_MRL_orients, ior, jor, kor);
00818        }
00819        else if( debug > 3 )
00820          fprintf(stderr,"-d dicom: bad orient vectors qq = %d, xn,yn = %f,%f\n",
00821                  qq, xn, yn);
00822      }
00823 
00824    } else if( epos[E_PATIENT_ORIENTATION] != NULL ){
00825               /* symbolic orientation of image [not so useful, or common] */
00826      ddd = strstr(epos[E_PATIENT_ORIENTATION],"//") ;
00827      if( ddd != NULL ){
00828        char xc='\0' , yc='\0' ;
00829        sscanf(ddd+2,"%c\\%c",&xc,&yc) ;   /* e.g., "L\P" */
00830        switch( toupper(xc) ){
00831          case 'L': DI_MRL_orients[0]='L'; DI_MRL_orients[1]='R'; ior=-1; break;
00832          case 'R': DI_MRL_orients[0]='R'; DI_MRL_orients[1]='L'; ior= 1; break;
00833          case 'P': DI_MRL_orients[0]='P'; DI_MRL_orients[1]='A'; ior=-2; break;
00834          case 'A': DI_MRL_orients[0]='A'; DI_MRL_orients[1]='P'; ior= 2; break;
00835          case 'F': DI_MRL_orients[0]='I'; DI_MRL_orients[1]='S'; ior= 3; break;
00836                    /* F=foot */
00837          case 'H': DI_MRL_orients[0]='S'; DI_MRL_orients[1]='I'; ior=-3; break;
00838                    /* H=head */
00839          default:  DI_MRL_orients[0]='\0';DI_MRL_orients[1]='\0'; ior= 0; break;
00840        }
00841        switch( toupper(yc) ){
00842          case 'L': DI_MRL_orients[2]='L'; DI_MRL_orients[3]='R'; jor=-1; break;
00843          case 'R': DI_MRL_orients[2]='R'; DI_MRL_orients[3]='L'; jor= 1; break;
00844          case 'P': DI_MRL_orients[2]='P'; DI_MRL_orients[3]='A'; jor=-2; break;
00845          case 'A': DI_MRL_orients[2]='A'; DI_MRL_orients[3]='P'; jor= 2; break;
00846          case 'F': DI_MRL_orients[2]='I'; DI_MRL_orients[3]='S'; jor= 3; break;
00847          case 'H': DI_MRL_orients[2]='S'; DI_MRL_orients[3]='I'; jor=-3; break;
00848          default:  DI_MRL_orients[2]='\0';DI_MRL_orients[3]='\0'; jor= 0; break;
00849        }
00850        DI_MRL_orients[6]='\0' ;      /* terminate orientation string */
00851        kor = 6 - abs(ior)-abs(jor) ;   /* which spatial direction is z-axis */
00852        have_orients = (ior != 0 && jor != 0) ;
00853 
00854        if( debug > 3 )
00855          fprintf(stderr,"-d dicom: DI_MRL_orients P = %s, [ijk]or = %d,%d,%d\n",
00856                  DI_MRL_orients, ior, jor, kor);
00857      }
00858 
00859    }  /* end of 2D image orientation */
00860 
00861    if( !have_orients ){
00862       fprintf(stderr,"** MRD: failed to determine dicom image orientation\n");
00863       free(ppp);  free(im);
00864       if( data ){ free(*data); *data = NULL; }
00865       RETURN(NULL);
00866    }
00867 
00868    /* skip mosaic use for finishing orientation string */
00869 
00870    /** use image position vector to set offsets,
00871        and (2cd time in) the z-axis orientation **/
00872 
00873    if( epos[E_IMAGE_POSITION] == NULL ){
00874       fprintf(stderr,"** MRD: missing image position\n");
00875       free(ppp);  free(im);
00876       if( data ){ free(*data); *data = NULL; }
00877       RETURN(NULL);
00878    }
00879 
00880    ddd = strstr(epos[E_IMAGE_POSITION],"//") ;
00881    if( ddd == NULL ){
00882       fprintf(stderr,"** MRD: missing IMAGE_POSITION\n");
00883       free(ppp);  free(im);
00884       if( data ){ free(*data); *data = NULL; }
00885       RETURN(NULL);
00886    }
00887 
00888    {   /* old mri_read_dicom: line 982 */
00889        float xyz[3] ; int qq ;
00890        qq = sscanf(ddd+2,"%f\\%f\\%f",xyz,xyz+1,xyz+2) ;
00891        if( qq != 3 ){
00892          fprintf(stderr,"** MRD: failed to read IMAGE_POSITION\n");
00893          free(ppp); free(im);
00894          if( data ){ free(*data); *data = NULL; }
00895          RETURN(NULL);
00896        }
00897        /* now use [ijk]or from above */
00898        im->xo = xyz[abs(ior)-1];
00899        im->yo = xyz[abs(jor)-1];
00900        im->zo = xyz[abs(kor)-1];
00901 
00902        if( debug > 3 )
00903           fprintf(stderr,"-d dicom: read RAI image position %f, %f, %f\n"
00904                          "          and dset image position %f, %f, %f\n",
00905                   xyz[0], xyz[1], xyz[2], im->xo, im->yo, im->zo );
00906 
00907        /* fill the orients string */
00908        {
00909          static float zoff ;      /* saved from nzoff=0 case */
00910          float zz = xyz[kor-1] ;  /* kor from orients above */
00911 
00912          if( nzoff == 0 ){  /* 1st DICOM image */
00913 
00914            zoff = zz ;      /* save this for 2nd image calculation */
00915 
00916            /* 01 Nov 2002: in mosaic case, may have set this already */
00917 
00918            if( DI_MRL_orients[4] == '\0' ){
00919              switch( kor ){   /* may be changed on second image */
00920                case  1: DI_MRL_orients[4] = 'L'; DI_MRL_orients[5] = 'R'; break;
00921                case  2: DI_MRL_orients[4] = 'P'; DI_MRL_orients[5] = 'A'; break;
00922                case  3: DI_MRL_orients[4] = 'I'; DI_MRL_orients[5] = 'S'; break;
00923                default: DI_MRL_orients[4] ='\0'; DI_MRL_orients[5] ='\0'; break;
00924              }
00925              DI_MRL_orients[6] = '\0' ;
00926            }
00927 
00928            /* Save x,y offsets of this 1st slice */
00929 
00930            qq = abs(ior) ;
00931            DI_MRL_xoff = xyz[qq-1] ; use_DI_MRL_xoff = 1 ;
00932            if( ior > 0 ) DI_MRL_xoff = -DI_MRL_xoff ;
00933 
00934            qq = abs(jor) ;
00935            DI_MRL_yoff = xyz[qq-1] ; use_DI_MRL_yoff = 1 ;
00936            if( jor > 0 ) DI_MRL_yoff = -DI_MRL_yoff ;
00937 
00938          } else if( nzoff == 1 && !use_DI_MRL_zoff ){  /* 2nd DICOM image */
00939 
00940            float qoff = zz - zoff ;    /* vive la difference */
00941            if( qoff < 0 ) kor = -kor ; /* kor determines z-axis orientation */
00942 #if 0
00943 fprintf(stderr,"  nzoff=1 kor=%d qoff=%f\n",kor,qoff) ;
00944 #endif
00945            switch( kor ){
00946              case  1: DI_MRL_orients[4] = 'R'; DI_MRL_orients[5] = 'L'; break;
00947              case -1: DI_MRL_orients[4] = 'L'; DI_MRL_orients[5] = 'R'; break;
00948              case  2: DI_MRL_orients[4] = 'A'; DI_MRL_orients[5] = 'P'; break;
00949              case -2: DI_MRL_orients[4] = 'P'; DI_MRL_orients[5] = 'A'; break;
00950              case  3: DI_MRL_orients[4] = 'I'; DI_MRL_orients[5] = 'S'; break;
00951              case -3: DI_MRL_orients[4] = 'S'; DI_MRL_orients[5] = 'I'; break;
00952              default: DI_MRL_orients[4] ='\0'; DI_MRL_orients[5] ='\0'; break;
00953            }
00954            DI_MRL_orients[6] = '\0' ;
00955 
00956            /* save spatial offset of first slice              */
00957            /* [this needs to be positive in the direction of] */
00958            /* [the -z axis, so may need to change its sign  ] */
00959 
00960            DI_MRL_zoff = zoff ; use_DI_MRL_zoff = 1 ;
00961            if( kor > 0 ) DI_MRL_zoff = -DI_MRL_zoff ;
00962          }
00963          nzoff++ ;  /* 3rd and later images don't count for z-orientation */
00964        }
00965    }  /* end of using image position */
00966 
00967    /** 23 Dec 2002:
00968        use slice location value to set z-offset,
00969        and (2cd time in) the z-axis orientation
00970        -- only try this if image position vector (code above) isn't present
00971           AND if we don't have a mosaic image (which already did this stuff)
00972        -- shouldn't be necessary, since slice location is deprecated        **/
00973 
00974    /* if this is here, use it for additional accuracy
00975       (it must basically match current zo) */
00976 
00977    { /* just use one warning counter, since these should be unique */
00978      static int nwarn = 0;
00979 
00980      if( ( epos[E_SLICE_LOCATION] == NULL) ||
00981          ( (ddd = strstr(epos[E_SLICE_LOCATION],"//")) == NULL ) )
00982      {
00983        if( nwarn == 0 && debug > 1 )
00984           fprintf(stderr,"-d dimon: missing SLICE_LOCATION, continuing...\n");
00985        nwarn++;
00986      } else {
00987         /* get and test the slice location */
00988         float zz ; int qq;
00989         qq = sscanf(ddd+2,"%f",&zz) ;
00990         if( qq != 1 ){
00991            if( !nwarn )fprintf(stderr,"** failed to extract SLICE_LOCATION\n");
00992            nwarn++;
00993         } else {
00994            /* it seems we have to add our own sign to the slice location */
00995            if( zz * im->zo < 0.0 ){
00996               if( (nwarn == 0) && (debug > 3) )
00997                 fprintf(stderr,"-d image and slice loc diff in sign: %f, %f\n",
00998                         im->zo, zz);
00999               nwarn++;
01000               zz = -zz;
01001            }
01002   
01003            if( fabs(zz - im->zo) > 0.1 ){
01004               fprintf(stderr,
01005                       "** MRD: IMAGE_LOCATION and SLICE_LOCATION disagree!\n"
01006                       "   z coord from IL = %f, from SL = %f\n", im->zo,zz);
01007               free(ppp); free(im);
01008               if( data ){ free(*data); *data = NULL; }
01009               RETURN(NULL);
01010            }
01011 
01012            if( debug > 3 )
01013               fprintf(stderr,"-d dicom: using slice location %f (zo = %f)\n",
01014                       zz, im->zo );
01015            im->zo = zz;
01016         }
01017      }
01018    }
01019 
01020    free(ppp);  /* free the ASCII header */
01021 
01022    RETURN( im );
01023 }

00111 {
00112    long val = 1;
00113 
00114    LITTLE_ENDIAN_ARCHITECTURE = (*(char *)&val == 1);
00115 }

00086 {
00087    char * dp0 = (char *)data;
00088    char * dp1 = dp0 + 1;
00089    char * dp2 = dp0 + 2;
00090    char * dp3 = dp0 + 3;
00091    int    c;
00092 
00093    for( c = 0; c < nvals; c++, dp0 += 4, dp1 += 4, dp2 += 4, dp3 += 4 )
00094    {
00095       *dp0 ^= *dp3; *dp3 ^= *dp0; *dp0 ^= *dp3;
00096       *dp1 ^= *dp2; *dp2 ^= *dp1; *dp1 ^= *dp2;
00097    }
00098 }

00072 {
00073    char * dp0 = (char *)data;
00074    char * dp1 = dp0 + 1;
00075    int    c;
00076 
00077    for( c = 0; c < nvals; c++, dp0 += 2, dp1 += 2 )
00078    {
00079       *dp0 ^= *dp1;
00080       *dp1 ^= *dp0;
00081       *dp0 ^= *dp1;
00082    }
00083 }