#include "mrilib.h"
#include "vecmat.h"
/*-----------------------------------------------------------------------------------*/
#undef DEBUG_ON
#define NWMAX 2 /* max number of warning message of each type to print */
#define NMOMAX 256 /* never seen one this big!!! */
#define EPSILON 0.0001
#define ALMOST(a,b) \
( fabs(a - b) < EPSILON)
typedef struct {
int good ; /* data in here is good? */
int have_data[3] ; /* do we have slices 0 and 1 in *
* each dimension to determine z-spacing? *
* added 25 Feb 2003 KRH */
int mosaic_num ; /* how many slices in 1 'image' */
float slice_xyz[NMOMAX][3] ; /* Sag, Cor, Tra coordinates */
} Siemens_extra_info ;
typedef struct {
THD_fvec3 xvec, yvec; /* Image Orientation fields */
THD_fvec3 dfpos1; /* image origin for first two slices*/
THD_fvec3 dfpos2;
THD_fvec3 del; /* voxel dimensions */
int mosaic; /* data is mosaic */
int mos_ix, mos_nx, mos_ny, mos_nslice; /* mosaic properties */
int nx, ny; /* overall mosaic dimensions */
float Tr_dicom[4][4]; /* transformation matrix */
float slice_xyz[2][3]; /* coordinates for 1st and last slices */
int mos_sliceinfo; /* flag for existence of coordinate info */
} oblique_info;
oblique_info obl_info;
/* mod -16 May 2007 */
/* compute Tr transformation matrix for oblique data */
static float *ComputeObliquity(oblique_info *obl_info);
static void Clear_obl_info(oblique_info *obl_info);
static void Fill_obl_info(oblique_info *obl_info, char **epos, Siemens_extra_info *siem);
void mri_read_dicom_reset_obliquity();
void mri_read_dicom_get_obliquity(float *);
static char * extract_bytes_from_file( FILE *fp, off_t start, size_t len, int strize ) ;
static void get_siemens_extra_info( char *str , Siemens_extra_info *mi ) ;
static float get_dz( char **epos);
static int CheckObliquity(float xc1, float xc2, float xc3, float yc1, float yc2, float yc3);
static int obl_info_set = 0;
/*-----------------------------------------------------------------------------------*/
/* Save the Siemens extra info string in case the caller wants to get it. */
static char *str_sexinfo=NULL ;
char * mri_dicom_sexinfo(void){ return str_sexinfo; } /* 23 Dec 2002 */
/*-----------------------------------------------------------------------------------*/
static int LITTLE_ENDIAN_ARCHITECTURE = -1 ;
static void RWC_set_endianosity(void)
{
union { unsigned char bb[2] ;
short ss ; } fred ;
if( LITTLE_ENDIAN_ARCHITECTURE < 0 ){
fred.bb[0] = 1 ; fred.bb[1] = 0 ;
LITTLE_ENDIAN_ARCHITECTURE = (fred.ss == 1) ;
}
}
/*-----------------------------------------------------------------------------------*/
static char *elist[] = {
"0018 0050" , /* Slice thickness */
"0018 0080" , /* Repetition time */
"0018 0088" , /* Spacing between slices */
"0018 1149" , /* Field of view */
"0020 0020" , /* Patient orientation */
"0020 0032" , /* Image position (patient) */
"0020 0037" , /* Image orientation (patient) */
"0020 1041" , /* Slice location */
"0028 0002" , /* Samples per pixel */
"0028 0008" , /* Number of frames */
"0028 0010" , /* Rows */
"0028 0011" , /* Columns */
"0028 0030" , /* Pixel spacing */
"0028 0100" , /* Bits allocated */
"0028 0101" , /* Bits stored */
"0028 1052" , /* Rescale intercept */
"0028 1053" , /* Rescale slope */
"0028 1054" , /* Rescale type */
"0028 0004" , /* Photometric interpretation */
"0028 0103" , /* Pixel representation */
"0028 0102" , /* High bit */
"0028 1050" , /* Window center */
"0028 1051" , /* Window width */
"0008 0008" , /* ID Image type */
"0008 0070" , /* ID Manufacturer */
"0018 1310" , /* Acquisition Matrix */
"0029 1010" , /* Siemens addendum #1 */
"0029 1020" , /* Siemens addendum #2 */
"0002 0010" , /* Transfer Syntax [RWC - 05 Jul 2006] */
NULL } ;
#define NUM_ELIST (sizeof(elist)/sizeof(char *)-1)
#define E_SLICE_THICKNESS 0
#define E_REPETITION_TIME 1
#define E_SLICE_SPACING 2
#define E_FIELD_OF_VIEW 3
#define E_PATIENT_ORIENTATION 4
#define E_IMAGE_POSITION 5
#define E_IMAGE_ORIENTATION 6
#define E_SLICE_LOCATION 7
#define E_SAMPLES_PER_PIXEL 8
#define E_NUMBER_OF_FRAMES 9
#define E_ROWS 10
#define E_COLUMNS 11
#define E_PIXEL_SPACING 12
#define E_BITS_ALLOCATED 13
#define E_BITS_STORED 14
#define E_RESCALE_INTERCEPT 15
#define E_RESCALE_SLOPE 16
#define E_RESCALE_TYPE 17
#define E_PHOTOMETRIC_INTERPRETATION 18
#define E_PIXEL_REPRESENTATION 19
#define E_HIGH_BIT 20
#define E_WINDOW_CENTER 21
#define E_WINDOW_WIDTH 22
#define E_ID_IMAGE_TYPE 23 /* 28 Oct 2002: for Siemens mosaic */
#define E_ID_MANUFACTURER 24
#define E_ACQ_MATRIX 25
#define E_SIEMENS_1 26 /* 31 Oct 2002 */
#define E_SIEMENS_2 27
#define E_TRANSFER_SYNTAX 28 /* 05 Jul 2006 */
/*---------------------------------------------------------------------------*/
/* global set via 'to3d -assume_dicom_mosaic' 13 Mar 2006 [rickr] */
int assume_dicom_mosaic = 0; /* (case of 1 is equivalent to Rich's change) */
/*---------------------------------------------------------------------------*/
/*! Read image(s) from a DICOM file, if possible.
-----------------------------------------------------------------------------*/
MRI_IMARR * mri_read_dicom( char *fname )
{
char *ppp , *ddd ;
off_t poff ;
int plen ;
char *epos[NUM_ELIST] ;
int ii,jj , ee , bpp , datum ;
int nx,ny,nz , swap , shift=0 ;
float dx,dy,dz,dt ;
MRI_IMARR *imar ;
MRI_IMAGE *im ;
void *iar ;
FILE *fp ;
int have_orients=0 ;
int ior,jor,kor ;
static int nzoff=0 ; /* for determining z-axis orientation/offset from multiple calls */
int mosaic=0 , mos_nx,mos_ny , mos_ix,mos_iy,mos_nz ; /* 28 Oct 2002 */
Siemens_extra_info sexinfo ; /* 31 Oct 2002 */
#if 0
short sbot,stop ;
float xcen,ycen,zcen ;
int use_xycen=0 ;
#endif
float dxx,dyy,dzz ;
char *eee ;
float rescale_slope=0.0 , rescale_inter=0.0 ; /* 23 Dec 2002 */
float window_center=0.0 , window_width =0.0 ;
char *sexi_start, *sexi_start2; /* KRH 25 Jul 2003 */
char *sexi_end;
int ts_endian = 1 ; /* 05 Jul 2006 - transfer syntax endian-ness
- 1 = little endian, 0 = big - RWCox */
int un16 = 0 ; /* 05 Jul 2006 - is it 16 bit unsigned data? */
int ov16 = 0 ; /* - did 16 bit overflow occur? */
static int obliqueflag = 0;
float xc1=0.0,xc2=0.0,xc3=0.0 , yc1=0.0,yc2=0.0,yc3=0.0 ;
float xn,yn ; int qq ;
ENTRY("mri_read_dicom") ;
if( str_sexinfo != NULL ){ free(str_sexinfo); str_sexinfo=NULL; }
if( !mri_possibly_dicom(fname) ){ /* 07 May 2003 */
if( AFNI_yesenv("AFNI_DICOM_VERBOSE") )
ERROR_message("file %s is not possibly DICOM",fname) ;
RETURN(NULL) ;
}
/* extract header info from file into a string
- cf. mri_dicom_hdr.[ch]
- run 'dicom_hdr -noname fname' to see the string format */
sexinfo.mosaic_num = 1; /* initialize to non-mosaic */
mri_dicom_nohex(1) ; /* don't print ints in hex mode */
mri_dicom_noname(1) ; /* don't print names, just tags */
ppp = mri_dicom_header( fname ) ; /* print header to malloc()-ed string */
if( ppp == NULL ){ /* didn't work; not a DICOM file? */
if( AFNI_yesenv("AFNI_DICOM_VERBOSE") )
ERROR_message("file %s is not interpretable as DICOM",fname) ;
RETURN(NULL) ;
}
/* find positions in header of elements we care about */
for( ee=0 ; ee < NUM_ELIST ; ee++ )
epos[ee] = strstr(ppp,elist[ee]) ;
/* Determine if we need to swap bytes after reading image data */
/* DICOM files are stored in LSB first (little endian) mode */
/* 05 Jul 2006 - not necessarily: check transfer syntax - RWC */
RWC_set_endianosity() ;
/* Transfer Syntax possibilities for the image data are:
"1.2.840.10008.1.2" = implicit little endian (default)
"1.2.840.10008.1.2.1" = explicit little endian
"1.2.840.10008.1.2.2" = explicit big endian
Various other possibilities (not supported herein) include:
"1.2.840.10008.1.2.4.70" = lossless JPEG
"1.2.840.10008.1.2.4.57" = another lossless JPEG
"1.2.840.10008.1.2.4.50" = 8-bit lossy JPEG
"1.2.840.10008.1.2.4.51" = 12-bit lossy JPEG
"1.2.840.10008.1.2.5" = RLE compression
"1.2.840.10008.1.2.4.80" = lossless JPEG-LS
"1.2.840.10008.1.2.4.81" = near-lossless JPEG-LS
"1.2.840.10008.1.2.4.90" = lossless JPEG-2000
"1.2.840.10008.1.2.4.91" = lossy JPEG-2000
"1.2.840.10008.1.2.4.100" = MPEG-2
"1.2.840.10008.1.2.1.99" = 'deflate' compression
Plus, 'private' transfer syntaxes are legal. 05 Jul 2006 - RWCox */
if( epos[E_TRANSFER_SYNTAX] != NULL ){
ddd = strstr(epos[E_TRANSFER_SYNTAX],"//") ;
if( ddd != NULL ){
char ts[256]="\0" ;
sscanf(ddd+2,"%254s",ts) ;
ts_endian = -1 ;
if( strlen(ts) >= 17 && strncmp(ts,"1.2.840.10008.1.2",17)==0 ){
if( strcmp(ts,"1.2.840.10008.1.2.2") == 0 ) /* big endian */
ts_endian = 0 ; /* big endian */
else if( strcmp(ts,"1.2.840.10008.1.2.1") == 0 ||
strcmp(ts,"1.2.840.10008.1.2" ) == 0 ) /* little endian */
ts_endian = 1 ;
}
if( ts_endian < 0 ){
static int nwarn=0 ;
if( nwarn < NWMAX )
WARNING_message("DICOM file %s: unsupported Transfer Syntax '%s'",fname,ts) ;
if( nwarn == NWMAX )
WARNING_message("DICOM: no more Transfer Syntax messages will be printed") ;
nwarn++ ;
}
}
}
if( ts_endian < 0 ) ts_endian = 1 ;
swap = (ts_endian != LITTLE_ENDIAN_ARCHITECTURE) ;
/* find out where the pixel array is in the file */
mri_dicom_pxlarr( &poff , &plen ) ;
if( plen <= 0 || poff == 0 ){
ERROR_message("DICOM file %s: ILLEGAL plen=%d poff=%u",
fname , plen , (unsigned int)poff ) ;
free(ppp); RETURN(NULL);
}
/* check if file is actually this big (maybe it was truncated) */
{ unsigned int psiz , fsiz ;
fsiz = THD_filesize( fname ) ;
psiz = (unsigned int)(poff) + (unsigned int)(plen) ;
if( fsiz < psiz ){
ERROR_message("DICOM file %s: plen=%d poff=%u filesize=%u",
fname , plen , (unsigned int)poff , fsiz ) ;
free(ppp) ; RETURN(NULL);
}
}
/* see if the header has the elements we absolutely need */
if( epos[E_ROWS] == NULL ||
epos[E_COLUMNS] == NULL ||
epos[E_BITS_ALLOCATED] == NULL ){
ERROR_message("DICOM file %s: missing mandatory fields",fname) ;
free(ppp) ; RETURN(NULL);
}
/* check if we have 1 sample per pixel (can't deal with 3 or 4 now) */
if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
if( ii != 1 ){
ERROR_message("DICOM file %s: %d samples per pixel -- too much!",
fname,ii) ;
free(ppp) ; RETURN(NULL);
}
}
/* check if photometric interpretation is MONOCHROME (don't like PALETTE) */
if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
if( ddd == NULL ){
ERROR_message("DICOM file %s: not MONOCHROME!",fname) ;
free(ppp) ; RETURN(NULL);
}
}
/* check if we have 8, 16, or 32 bits per pixel */
ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
if( ddd == NULL ){ free(ppp); RETURN(NULL); }
bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
switch( bpp ){
default:
ERROR_message(
"DICOM file %s: unsupported Bits Allocated = %d", fname , bpp ) ;
free(ppp); RETURN(NULL); /* bad value */
case 8: datum = MRI_byte ; break ;
case 16: datum = MRI_short; break ;
case 32: datum = MRI_int ; break ; /* probably not present in DICOM? */
}
bpp /= 8 ; /* now bytes per pixel, instead of bits */
/*** Print some warnings if appropriate ***/
/* check if BITS_STORED and HIGH_BIT are aligned */
if( epos[E_BITS_STORED] != NULL && epos[E_HIGH_BIT] != NULL ){
int bs=0 , hb=0 ;
ddd = strstr(epos[E_BITS_STORED],"//") ; sscanf(ddd+2,"%d",&bs) ;
ddd = strstr(epos[E_HIGH_BIT],"//") ; sscanf(ddd+2,"%d",&hb) ;
if( bs != hb+1 ){
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,
"++ DICOM WARNING: file %s has Bits_Stored=%d and High_Bit=%d\n",
fname,bs,hb) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM WARNING: no more Bits_Stored messages will be printed\n") ;
nwarn++ ;
}
}
/* check if Rescale is ordered */
/* 23 Dec 2002: actually get the rescale params, if environment says to */
eee = getenv("AFNI_DICOM_RESCALE") ;
if( epos[E_RESCALE_INTERCEPT] != NULL && epos[E_RESCALE_SLOPE] != NULL ){
if( eee == NULL || toupper(*eee) != 'Y' ){
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,
"++ DICOM WARNING: file %s has Rescale tags; setenv AFNI_DICOM_RESCALE YES to enforce them\n",
fname ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM WARNING: no more Rescale tags messages will be printed\n") ;
nwarn++ ;
} else {
ddd = strstr(epos[E_RESCALE_INTERCEPT],"//") ; sscanf(ddd+2,"%f",&rescale_inter) ;
ddd = strstr(epos[E_RESCALE_SLOPE ],"//") ; sscanf(ddd+2,"%f",&rescale_slope) ;
}
}
/* check if Window is ordered */
/* 23 Dec 2002: actually get the window params, if environment says to */
eee = getenv("AFNI_DICOM_WINDOW") ;
if( epos[E_WINDOW_CENTER] != NULL && epos[E_WINDOW_WIDTH] != NULL ){
if( eee == NULL || toupper(*eee) != 'Y' ){
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,
"++ DICOM WARNING: file %s has Window tags; setenv AFNI_DICOM_WINDOW YES to enforce them\n",
fname ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM WARNING: no more Window tags messages will be printed\n") ;
nwarn++ ;
} else {
ddd = strstr(epos[E_WINDOW_CENTER],"//") ; sscanf(ddd+2,"%f",&window_center) ;
ddd = strstr(epos[E_WINDOW_WIDTH ],"//") ; sscanf(ddd+2,"%f",&window_width ) ;
}
}
/*** extract attributes of the image(s) to be read in ***/
/* get image nx & ny */
ddd = strstr(epos[E_ROWS],"//") ; /* 31 Oct 2002: */
if( ddd == NULL ){ free(ppp) ; RETURN(NULL); } /* Oops: ROWS is ny and */
ny = 0 ; sscanf(ddd+2,"%d",&ny) ; /* COLUMNS is nx! */
if( ny < 2 ){ free(ppp) ; RETURN(NULL); }
ddd = strstr(epos[E_COLUMNS],"//") ;
if( ddd == NULL ){ free(ppp) ; RETURN(NULL); }
nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
if( nx < 2 ){ free(ppp) ; RETURN(NULL); }
/* get number of slices */
nz = 0 ;
if( epos[E_NUMBER_OF_FRAMES] != NULL ){
ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
}
/* if didn't get nz above, make up a value */
if( nz == 0 ) nz = plen / (bpp*nx*ny) ; /* compute from image array size */
if( nz == 0 ){
ERROR_message("DICOM file %s: not enough data for 1 slice!",fname) ;
free(ppp) ; RETURN(NULL);
}
/*-- 28 Oct 2002: Check if this is a Siemens mosaic. --*/
/*-- 02 Dec 2002: Don't use Acquisition Matrix anymore;
instead, use the Siemens extra info
in epos[E_SIEMENS_2]. --*/
/*-- 24 Dec 2002: Extract str_sexinfo even if not a mosaic. --*/
if( epos[E_ID_MANUFACTURER] != NULL &&
strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
epos[E_SIEMENS_2] != NULL ){
int len=0,loc=0 , aa,bb ;
sscanf(epos[E_SIEMENS_2],"%x%x%d [%d" , &aa,&bb , &len,&loc ) ;
if( len > 0 && loc > 0 ){
fp = fopen( fname , "rb" ) ;
if( fp != NULL ){
str_sexinfo = extract_bytes_from_file( fp, (off_t)loc, (size_t)len, 1 ) ;
fclose(fp) ;
}
}
}
/*-- process str_sexinfo only if this is marked as a mosaic image --*/
/*-- preliminary processing of sexinfo EVEN IF NOT MARKED AS MOSAIC, --*/
/*-- SINCE PSYCHOTIC, UNMARKED MOSAICS ARE TURNING UP IN THE WILD! KRH, 11/6/05 --*/
/* 31 Oct 2002: extract extra Siemens info from str_sexinfo */
/* KRH 25 Jul 2003 if start and end markers are present for
* Siemens extra info, cut string down to those boundaries */
/* if assume_dicom_mosaic is not set, then require "MOSAIC" string */
/* 13 Mar 2006 [rickr] */
if( ( assume_dicom_mosaic ||
( epos[E_ID_IMAGE_TYPE] && strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC") ) )
&& str_sexinfo != NULL ){
sexi_start = strstr(str_sexinfo, "### ASCCONV BEGIN ###");
if(sexi_start != NULL) { /* search for end after start - drg,fredtam 23 Mar 2007 */
sexi_start2 = strstr(sexi_start+21, "### ASCCONV BEGIN ###");
sexi_end = strstr(sexi_start, "### ASCCONV END ###");
if (sexi_end != NULL) {
char *sexi_tmp;
int sexi_size;
if((sexi_start2!=NULL) && (sexi_start2<sexi_end)) {
sexi_start = sexi_start2;
}
sexi_size = sexi_end - sexi_start + 19 ;
sexi_tmp = AFMALL( char, sexi_size );
memcpy(sexi_tmp,sexi_start,sexi_size);
free(str_sexinfo);
str_sexinfo = sexi_tmp;
}
}
/* end KRH 25 Jul 2003 change */
sexinfo.good = 0 ; /* start by marking it as bad */
for(ii = 0; ii < 3; ii++) sexinfo.have_data[ii] = 0; /* 25 Feb 03 Initialize new member KRH */
get_siemens_extra_info( str_sexinfo , &sexinfo ) ;
if( sexinfo.good ){ /* if data is good */
if(( epos[E_ID_IMAGE_TYPE] != NULL &&
strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC") != NULL) ||
sexinfo.mosaic_num > 1 ) {
/* compute size of mosaic layout
as 1st integer whose square is >= # of images in mosaic */
for( mos_ix=1 ; mos_ix*mos_ix < sexinfo.mosaic_num ; mos_ix++ ) ; /* nada */
mos_iy = mos_ix ;
mos_nx = nx / mos_ix ; mos_ny = ny / mos_iy ; /* sub-image dimensions */
if( mos_ix*mos_nx == nx && /* Sub-images must fit nicely */
mos_iy*mos_ny == ny ){ /* into super-image layout. */
static int nwarn=0 ;
/* should be tagged as a 1 slice image thus far */
if( nz > 1 ){
fprintf(stderr,
"** DICOM ERROR: %dx%d Mosaic of %dx%d images in file %s, but also have nz=%d\n",
mos_ix,mos_iy,mos_nx,mos_ny,fname,nz) ;
free(ppp) ; RETURN(NULL) ;
}
/* mark as a mosaic */
mosaic = 1 ;
mos_nz = mos_ix * mos_iy ; /* number of slices in mosaic */
if( nwarn < NWMAX )
fprintf(stderr,"++ DICOM NOTICE: %dx%d Siemens Mosaic of %d %dx%d images in file %s\n",
mos_ix,mos_iy,sexinfo.mosaic_num,mos_nx,mos_ny,fname) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic messages will be printed\n") ;
nwarn++ ;
} /* end of if mosaic sizes are reasonable */
else { /* warn about bad mosaic sizes */
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,
"++ DICOM WARNING: bad Siemens Mosaic params: nx=%d ny=%d ix=%d iy=%d imx=%d imy=%d\n",
mos_nx,mos_ny , mos_ix,mos_iy , nx,ny ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic param messages will be printed\n");
nwarn++ ;
}
} /* end of if a Siemens mosaic */
} /* end of if sexinfo was good */
else { /* warn if sexinfo was bad */
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,"++ DICOM WARNING: indecipherable Siemens Mosaic info (%s) in file %s\n",
elist[E_SIEMENS_2] , fname ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic info messages will be printed\n");
nwarn++ ;
}
} /* end of if str_sexinfo exists */
/*-- try to get dx, dy, dz, dt --*/
dx = dy = dz = dt = 0.0 ;
/* dx,dy first */
if( epos[E_PIXEL_SPACING] != NULL ){
ddd = strstr(epos[E_PIXEL_SPACING],"//") ;
if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
if( dy == 0.0 && dx > 0.0 ) dy = dx ;
}
if( dx == 0.0 && epos[E_FIELD_OF_VIEW] != NULL ){
ddd = strstr(epos[E_FIELD_OF_VIEW],"//") ;
if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
if( dx > 0.0 ){
if( dy == 0.0 ) dy = dx ;
dx /= nx ; dy /= ny ;
}
}
/*-- 27 Nov 2002: fix stupid GE error,
where the slice spacing is really the slice gap --*/
dz = get_dz(epos);
/* get dt */
if( epos[E_REPETITION_TIME] != NULL ){
ddd = strstr(epos[E_REPETITION_TIME],"//") ;
if( ddd != NULL ){
sscanf( ddd+2 , "%f" , &dt ) ;
dt *= 0.001 ; /* ms to s */
}
}
/* check if we might have 16 bit unsigned data that fills all bits */
if( bpp == 2 ){
if( epos[E_PIXEL_REPRESENTATION] != NULL ){
ddd = strstr(epos[E_PIXEL_REPRESENTATION],"//") ;
if( ddd != NULL ){
ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
if( ii == 0 ) un16 = 1 ; /* unsigned */
}
}
if( un16 && epos[E_HIGH_BIT] != NULL ){
ddd = strstr(epos[E_HIGH_BIT],"//") ;
if( ddd != NULL ){
ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
if( ii < 15 ) un16 = 0 ; /* but less than 16 bits */
}
}
}
#if 0
if( bpp == 2 ){
if( epos[E_PIXEL_REPRESENTATION] != NULL ){
ddd = strstr(epos[E_PIXEL_REPRESENTATION],"//") ;
if( ddd != NULL ){
ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
if( ii == 1 ) shift = -1 ;
}
}
if( shift == 0 && epos[E_HIGH_BIT] != NULL ){
ddd = strstr(epos[E_HIGH_BIT],"//") ;
if( ddd != NULL ){
ii = 0 ; sscanf( ddd+2 , "%d" , &ii ) ;
if( ii == 15 ) shift = 1 ;
}
}
sbot = 32767 ; stop = -32767 ;
}
#endif
/** Finally! Read images from file. **/
fp = fopen( fname , "rb" ) ;
if( fp == NULL ){
ERROR_message("DICOM file %s: can't open!?",fname) ;
free(ppp) ; RETURN(NULL);
}
lseek( fileno(fp) , poff , SEEK_SET ) ;
INIT_IMARR(imar) ;
/* 28 Oct 2002: must allow for 2D mosaic mode */
if( !mosaic ){ /*-- 28 Oct 2002: old method, not a mosaic --*/
for( ii=0 ; ii < nz ; ii++ ){
im = mri_new( nx , ny , datum ) ; /* new MRI_IMAGE struct */
iar = mri_data_pointer( im ) ; /* data array in struct */
fread( iar , bpp , nx*ny , fp ) ; /* read data directly into it */
if( swap ){ /* swap bytes? */
switch( im->pixel_size ){
default: break ;
case 2: swap_twobytes ( im->nvox, iar ) ; break ; /* short */
case 4: swap_fourbytes( im->nvox, iar ) ; break ; /* int, float */
case 8: swap_fourbytes( 2*im->nvox, iar ) ; break ; /* complex */
}
im->was_swapped = 1 ;
}
#if 0
if( shift == 1 ){
switch( datum ){
case MRI_short:{
short * sar = (short *) iar ;
for( jj=0 ; jj < im->nvox ; jj++ ){
sbot = MIN( sar[jj] , sbot ) ;
stop = MAX( sar[jj] , stop ) ;
}
}
break ;
}
}
#endif
/* store auxiliary data in image struct */
if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
}
if( dt > 0.0 ) im->dt = dt ;
ADDTO_IMARR(imar,im) ;
}
} else { /*-- 28 Oct 2002: is a 2D mosaic --*/
char *dar , *iar ;
int last_ii=-1 , nvox , yy,xx,nxx ;
nvox = mos_nx*mos_ny*mos_nz ; /* total number of voxels */
dar = (char*)calloc(bpp,nvox) ; /* make space for super-image */
fread( dar , bpp , nvox , fp ) ; /* read data directly into it */
if( swap ){ /* swap bytes? */
switch( bpp ){
default: break ;
case 2: swap_twobytes ( nvox, dar ) ; break ; /* short */
case 4: swap_fourbytes( nvox, dar ) ; break ; /* int, float */
case 8: swap_fourbytes( 2*nvox, dar ) ; break ; /* complex */
}
}
/* load data from dar into images */
nxx = mos_nx * mos_ix ; /* # pixels per mosaic line */
for( yy=0 ; yy < mos_iy ; yy++ ){ /* loop over sub-images in mosaic */
for( xx=0 ; xx < mos_ix ; xx++ ){
im = mri_new( mos_nx , mos_ny , datum ) ;
iar = mri_data_pointer( im ) ; /* sub-image array */
/* copy data rows from dar into iar */
for( jj=0 ; jj < mos_ny ; jj++ ) /* loop over rows inside sub-image */
memcpy( iar + jj*mos_nx*bpp ,
dar + xx*mos_nx*bpp + (jj+yy*mos_ny)*nxx*bpp ,
mos_nx*bpp ) ;
if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
}
if( dt > 0.0 ) im->dt = dt ;
if( swap ) im->was_swapped = 1 ;
ADDTO_IMARR(imar,im) ;
} /* end of x sub-image loop */
} /* end of y sub-image loop */
free(dar) ; /* don't need no more; copied all data out of it now */
/* truncate zero images out of tail of mosaic */
if( sexinfo.good ){ /* the new way: use the mosaic count from Siemens extra info */
if( sexinfo.mosaic_num < IMARR_COUNT(imar) )
TRUNCATE_IMARR(imar,sexinfo.mosaic_num) ;
} else { /* the old way: find the last image with a nonzero value inside */
for( ii=mos_nz-1 ; ii >= 0 ; ii-- ){ /* loop backwards over images */
im = IMARR_SUBIM(imar,ii) ;
switch( im->kind ){
case MRI_short:{
short *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
if( jj < im->nvox ) last_ii = ii ;
}
break ;
case MRI_byte:{
byte *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
if( jj < im->nvox ) last_ii = ii ;
}
break ;
case MRI_int:{
int *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ ) if( ar[jj] != 0 ) break ;
if( jj < im->nvox ) last_ii = ii ;
}
break ;
}
if( last_ii >= 0 ) break ;
}
if( last_ii <= 0 ) last_ii = 1 ;
if( last_ii+1 < IMARR_COUNT(imar) ) TRUNCATE_IMARR(imar,last_ii+1) ;
} /* end of truncating off all zero images at end */
#if 0
fprintf(stderr,"\nmri_read_dicom Mosaic: mos_nx=%d mos_ny=%d mos_ix=%d mos_iy=%d slices=%d\n",
mos_nx,mos_ny,mos_ix,mos_iy,IMARR_COUNT(imar)) ;
MCHECK ;
#endif
} /* end of mosaic input mode */
fclose(fp) ; /* 10 Sep 2002: oopsie - forgot to close file */
/*-- 05 Jul 2006 - check for 16 bit overflow --*/
if( un16 ){
for( ov16=ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
short *sar = MRI_SHORT_PTR( IMARR_SUBIM(imar,ii) ) ;
for( jj=0 ; jj < im->nvox ; jj++ ) if( sar[jj] < 0 ) ov16++ ;
}
if( ov16 )
WARNING_message(
"DICOM file %s: unsigned 16-bit; AFNI stores signed: %d pixels < 0" ,
fname , ov16 ) ;
}
/*-- 23 Dec 2002: implement Rescale, if ordered --*/
if( rescale_slope > 0.0 ){
for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
im = IMARR_SUBIM(imar,ii) ;
switch( im->kind ){
case MRI_byte:{
byte *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
}
break ;
case MRI_short:{
short *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
}
break ;
case MRI_int:{
int *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
}
break ;
}
}
} /* end of Rescale */
/*-- 23 Dec 2002: implement Window, if ordered --*/
/* section C.11.2.1.2 (page 503) */
if( window_width >= 1.0 ){
float wbot,wtop,wfac ;
int ymax=0 ;
/* get output range */
ddd = strstr(epos[E_BITS_STORED],"//") ;
if( ddd != NULL ){
ymax = 0 ; sscanf(ddd+2,"%d",&ymax) ;
if( ymax > 0 ) ymax = (1 << ymax) - 1 ;
}
if( ymax <= 0 ){
switch( IMARR_SUBIM(imar,0)->kind ){
case MRI_byte: ymax = MRI_maxbyte ; break ;
case MRI_short: ymax = MRI_maxshort ; break ;
case MRI_int: ymax = MRI_maxint ; break ;
}
}
/** window_width == 1 is special **/
if( window_width == 1.0 ){ /** binary threshold case **/
wbot = window_center - 0.5 ; /* the threshold */
for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
im = IMARR_SUBIM(imar,ii) ;
switch( im->kind ){
case MRI_byte:{
byte *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
}
break ;
case MRI_short:{
short *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
}
break ;
case MRI_int:{
int *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0 : ymax ;
}
break ;
}
}
} else { /** linear windowing case **/
wbot = (window_center - 0.5) - 0.5*(window_width-1.0) ;
wtop = (window_center - 0.5) + 0.5*(window_width-1.0) ;
wfac = ymax / (window_width-1.0) ;
for( ii=0 ; ii < IMARR_COUNT(imar) ; ii++ ){
im = IMARR_SUBIM(imar,ii) ;
switch( im->kind ){
case MRI_byte:{
byte *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0
:(ar[jj] > wtop) ? ymax
: wfac*(ar[jj]-wbot)+0.499 ;
}
break ;
case MRI_short:{
short *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0
:(ar[jj] > wtop) ? ymax
: wfac*(ar[jj]-wbot)+0.499 ;
}
break ;
case MRI_int:{
int *ar = mri_data_pointer( im ) ;
for( jj=0 ; jj < im->nvox ; jj++ )
ar[jj] = (ar[jj] <= wbot) ? 0
:(ar[jj] > wtop) ? ymax
: wfac*(ar[jj]-wbot)+0.499 ;
}
break ;
}
}
}
} /* end of Window */
/*-- store some extra information in MRILIB globals, too? --*/
if( dt > 0.0 && MRILIB_tr <= 0.0 ) MRILIB_tr = dt ; /* TR */
/*-- try to get image orientation fields (also, set ior,jor,kor) --*/
if( epos[E_IMAGE_ORIENTATION] != NULL ){ /* direction cosines of image plane */
ddd = strstr(epos[E_IMAGE_ORIENTATION],"//") ;
if( ddd != NULL ){
qq = sscanf(ddd+2,"%f\\%f\\%f\\%f\\%f\\%f",&xc1,&xc2,&xc3,&yc1,&yc2,&yc3);
xn = sqrt( xc1*xc1 + xc2*xc2 + xc3*xc3 ) ; /* vector norms */
yn = sqrt( yc1*yc1 + yc2*yc2 + yc3*yc3 ) ;
if( qq == 6 && xn > 0.0 && yn > 0.0 ){ /* both vectors OK */
xc1 /= xn ; xc2 /= xn ; xc3 /= xn ; /* normalize vectors */
yc1 /= yn ; yc2 /= yn ; yc3 /= yn ;
if(!obl_info_set) {
obliqueflag = CheckObliquity(xc1, xc2, xc3, yc1, yc2, yc3);
if(obliqueflag) {
INFO_message("Data detected to be oblique");
/* can also check obliquity consistency here across dicom files*/
}
}
if( !use_MRILIB_xcos ){
MRILIB_xcos[0] = xc1 ; MRILIB_xcos[1] = xc2 ; /* save direction */
MRILIB_xcos[2] = xc3 ; use_MRILIB_xcos = 1 ; /* cosine vectors */
}
if( !use_MRILIB_ycos ){
MRILIB_ycos[0] = yc1 ; MRILIB_ycos[1] = yc2 ;
MRILIB_ycos[2] = yc3 ; use_MRILIB_ycos = 1 ;
}
/* x-axis orientation */
/* ior determines which spatial direction is x-axis */
/* and is the direction that has the biggest change */
dxx = fabs(xc1) ; ior = 1 ;
dyy = fabs(xc2) ; if( dyy > dxx ){ ior=2; dxx=dyy; }
dzz = fabs(xc3) ; if( dzz > dxx ){ ior=3; }
dxx = MRILIB_xcos[ior-1] ; if( dxx < 0. ) ior = -ior;
if( MRILIB_orients[0] == '\0' ){
switch( ior ){
case -1: MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; break;
case 1: MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; break;
case -2: MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; break;
case 2: MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; break;
case 3: MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; break;
case -3: MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; break;
default: MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; break;
}
}
/* y-axis orientation */
/* jor determines which spatial direction is y-axis */
/* and is the direction that has the biggest change */
dxx = fabs(yc1) ; jor = 1 ;
dyy = fabs(yc2) ; if( dyy > dxx ){ jor=2; dxx=dyy; }
dzz = fabs(yc3) ; if( dzz > dxx ){ jor=3; }
dyy = MRILIB_ycos[jor-1] ; if( dyy < 0. ) jor = -jor;
if( MRILIB_orients[2] == '\0' ){
switch( jor ){
case -1: MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; break;
case 1: MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; break;
case -2: MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; break;
case 2: MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; break;
case 3: MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; break;
case -3: MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; break;
default: MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; break;
}
}
MRILIB_orients[6] = '\0' ; /* terminate orientation string */
kor = 6 - abs(ior)-abs(jor) ; /* which spatial direction is z-axis */
/* where 1=L-R, 2=P-A, 3=I-S */
have_orients = 1 ;
#if 0
fprintf(stderr,"MRILIB_orients=%s (from IMAGE_ORIENTATION)\n",MRILIB_orients) ;
#endif
}
}
} else if( epos[E_PATIENT_ORIENTATION] != NULL ){ /* symbolic orientation of image */
/* [not so useful, or common] */
ddd = strstr(epos[E_PATIENT_ORIENTATION],"//") ;
if( ddd != NULL ){
char xc='\0' , yc='\0' ;
sscanf(ddd+2,"%c\\%c",&xc,&yc) ; /* e.g., "L\P" */
switch( toupper(xc) ){
case 'L': MRILIB_orients[0] = 'L'; MRILIB_orients[1] = 'R'; ior=-1; break;
case 'R': MRILIB_orients[0] = 'R'; MRILIB_orients[1] = 'L'; ior= 1; break;
case 'P': MRILIB_orients[0] = 'P'; MRILIB_orients[1] = 'A'; ior=-2; break;
case 'A': MRILIB_orients[0] = 'A'; MRILIB_orients[1] = 'P'; ior= 2; break;
case 'F': MRILIB_orients[0] = 'I'; MRILIB_orients[1] = 'S'; ior= 3; break; /* F = foot */
case 'H': MRILIB_orients[0] = 'S'; MRILIB_orients[1] = 'I'; ior=-3; break; /* H = head */
default: MRILIB_orients[0] ='\0'; MRILIB_orients[1] ='\0'; ior= 0; break;
}
switch( toupper(yc) ){
case 'L': MRILIB_orients[2] = 'L'; MRILIB_orients[3] = 'R'; jor=-1; break;
case 'R': MRILIB_orients[2] = 'R'; MRILIB_orients[3] = 'L'; jor= 1; break;
case 'P': MRILIB_orients[2] = 'P'; MRILIB_orients[3] = 'A'; jor=-2; break;
case 'A': MRILIB_orients[2] = 'A'; MRILIB_orients[3] = 'P'; jor= 2; break;
case 'F': MRILIB_orients[2] = 'I'; MRILIB_orients[3] = 'S'; jor= 3; break;
case 'H': MRILIB_orients[2] = 'S'; MRILIB_orients[3] = 'I'; jor=-3; break;
default: MRILIB_orients[2] ='\0'; MRILIB_orients[3] ='\0'; jor= 0; break;
}
MRILIB_orients[6] = '\0' ; /* terminate orientation string */
kor = 6 - abs(ior)-abs(jor) ; /* which spatial direction is z-axis */
have_orients = (ior != 0 && jor != 0) ;
}
} /* end of 2D image orientation */
/*-- try to get image offset (position), if have orientation from above --*/
if( nzoff == 0 && have_orients && mosaic && sexinfo.good ){ /* 01 Nov 2002: use Siemens mosaic info */
int qq ;
float z0, z1 ;
/* 25 Feb 2003 changing error checking for mosaics missing one or more *
* dimension of slice coordinates KRH */
if (sexinfo.have_data[kor-1]) {
z0 = sexinfo.slice_xyz[0][kor-1] ; /* kor from orients above */
z1 = sexinfo.slice_xyz[1][kor-1] ; /* z offsets of 1st 2 slices */
} else { /* warn if sexinfo was bad */
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,"++ DICOM WARNING: Unusable coord. in Siemens Mosaic info (%s) in file %s\n",
elist[E_SIEMENS_2] , fname ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM NOTICE: no more Siemens Mosaic info messages will be printed\n");
nwarn++ ;
}
#if 0
/* Save x,y center of this 1st slice */
xcen = sexinfo.slice_xyz[0][abs(ior)-1] ;
ycen = sexinfo.slice_xyz[0][abs(jor)-1] ; use_xycen = 1 ;
#endif
/* finish z orientation now */
if( z1-z0 < 0.0 ) kor = -kor ; /* reversed orientation */
if( MRILIB_orients[4] == '\0' ){
switch( kor ){
case 1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
case 2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
}
}
MRILIB_orients[6] = '\0' ;
#if 0
fprintf(stderr,"z0=%f z1=%f kor=%d MRILIB_orients=%s\n",z0,z1,kor,MRILIB_orients) ;
#endif
MRILIB_zoff = z0 ; use_MRILIB_zoff = 1 ;
if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
}
/** use image position vector to set offsets,
and (2cd time in) the z-axis orientation **/
if( nzoff < 2 && epos[E_IMAGE_POSITION] != NULL && have_orients ){
ddd = strstr(epos[E_IMAGE_POSITION],"//") ;
if( ddd != NULL ){
float xyz[3] ; int qq ;
qq = sscanf(ddd+2,"%f\\%f\\%f",xyz,xyz+1,xyz+2) ;
if( qq == 3 ){
static float zoff ; /* saved from nzoff=0 case */
float zz = xyz[kor-1] ; /* kor from orients above */
#if 0
fprintf(stderr,"IMAGE_POSITION=%f %f %f kor=%d\n",xyz[0],xyz[1],xyz[2],kor) ;
#endif
if( nzoff == 0 ){ /* 1st DICOM image */
zoff = zz ; /* save this for 2nd image calculation */
/* 01 Nov 2002: in mosaic case, may have set this already */
if( MRILIB_orients[4] == '\0' ){
switch( kor ){ /* may be changed on second image */
case 1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
case 2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
}
MRILIB_orients[6] = '\0' ;
}
/* Save x,y offsets of this 1st slice */
qq = abs(ior) ;
MRILIB_xoff = xyz[qq-1] ; use_MRILIB_xoff = 1 ;
if( ior > 0 ) MRILIB_xoff = -MRILIB_xoff ;
qq = abs(jor) ;
MRILIB_yoff = xyz[qq-1] ; use_MRILIB_yoff = 1 ;
if( jor > 0 ) MRILIB_yoff = -MRILIB_yoff ;
/* 01 Nov 2002: adjust x,y offsets for mosaic */
if( mosaic ){
if( MRILIB_xoff < 0.0 ) MRILIB_xoff += 0.5*dx*mos_nx*(mos_ix-1) ;
else MRILIB_xoff -= 0.5*dx*mos_nx*(mos_ix-1) ;
if( MRILIB_yoff < 0.0 ) MRILIB_yoff += 0.5*dy*mos_ny*(mos_iy-1) ;
else MRILIB_yoff -= 0.5*dy*mos_ny*(mos_iy-1) ;
}
} else if( nzoff == 1 && !use_MRILIB_zoff ){ /* 2nd DICOM image */
float qoff = zz - zoff ; /* vive la difference */
if( qoff < 0 ) kor = -kor ; /* kor determines z-axis orientation */
#if 0
fprintf(stderr," nzoff=1 kor=%d qoff=%f\n",kor,qoff) ;
#endif
switch( kor ){
case 1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
case 2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
}
MRILIB_orients[6] = '\0' ;
/* save spatial offset of first slice */
/* [this needs to be positive in the direction of] */
/* [the -z axis, so may need to change its sign ] */
MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
}
nzoff++ ; /* 3rd and later images don't count for z-orientation */
}
}
} /* end of using image position */
/** 23 Dec 2002:
use slice location value to set z-offset,
and (2cd time in) the z-axis orientation
-- only try this if image position vector (code above) isn't present
AND if we don't have a mosaic image (which already did this stuff)
-- shouldn't be necessary, since slice location is deprecated **/
else if( nzoff < 2 && epos[E_SLICE_LOCATION] != NULL && have_orients && !mosaic ){
ddd = strstr(epos[E_SLICE_LOCATION],"//") ;
if( ddd != NULL ){
float zz ; int qq ;
qq = sscanf(ddd+2,"%f",&zz) ;
if( qq == 1 ){
static float zoff ; /* saved from nzoff=0 case */
#if 0
fprintf(stderr,"SLICE_LOCATION = %f\n",zz) ;
#endif
if( nzoff == 0 ){ /* 1st DICOM image */
zoff = zz ; /* save this for 2nd image calculation */
if( MRILIB_orients[4] == '\0' ){
switch( kor ){ /* may be changed on second image */
case 1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
case 2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
}
MRILIB_orients[6] = '\0' ;
}
} else if( nzoff == 1 && !use_MRILIB_zoff ){ /* 2nd DICOM image */
float qoff = zz - zoff ; /* vive la difference */
if( qoff < 0 ) kor = -kor ; /* kor determines z-axis orientation */
switch( kor ){
case 1: MRILIB_orients[4] = 'R'; MRILIB_orients[5] = 'L'; break;
case -1: MRILIB_orients[4] = 'L'; MRILIB_orients[5] = 'R'; break;
case 2: MRILIB_orients[4] = 'A'; MRILIB_orients[5] = 'P'; break;
case -2: MRILIB_orients[4] = 'P'; MRILIB_orients[5] = 'A'; break;
case 3: MRILIB_orients[4] = 'I'; MRILIB_orients[5] = 'S'; break;
case -3: MRILIB_orients[4] = 'S'; MRILIB_orients[5] = 'I'; break;
default: MRILIB_orients[4] ='\0'; MRILIB_orients[5] ='\0'; break;
}
MRILIB_orients[6] = '\0' ;
/* save spatial offset of first slice */
/* [this needs to be positive in the direction of] */
/* [the -z axis, so may need to change its sign ] */
MRILIB_zoff = zoff ; use_MRILIB_zoff = 1 ;
if( kor > 0 ) MRILIB_zoff = -MRILIB_zoff ;
}
nzoff++ ; /* 3rd and later images don't count for z-orientation */
}
}
} /* end of using slice location */
/* perhaps shift data shorts */
#if 0
if( shift == 1 ){
switch( datum ){
case MRI_short:{
if( sbot < 0 ){
unsigned short *sar ; int nvox = IMARR_SUBIM(imar,0)->nvox ;
for( ii=0 ; ii < nz ; ii++ ){
sar = mri_data_pointer( IMARR_SUBIM(imar,ii) ) ;
for( jj=0 ; jj < nvox ; jj++ ) sar[jj] >>= 1 ;
}
}
}
break ;
}
}
#endif
if(obl_info_set<2)
Fill_obl_info(&obl_info, epos, &sexinfo);
free(ppp); RETURN( imar );
}
/*---------- compute slice thickness from DICOM header ----------*/
static float get_dz( char **epos)
{
int stupid_ge_fix , no_stupidity ;
float sp=0.0 , th=0.0, dz = 0.0 ;
static int nwarn=0 ;
char *eee, *ddd ;
eee = getenv("AFNI_SLICE_SPACING_IS_GAP") ;
stupid_ge_fix = (eee != NULL && (*eee=='Y' || *eee=='y') ) ;
no_stupidity = (eee != NULL && (*eee=='N' || *eee=='n') ) ; /* 03 Mar 2003 */
if( epos[E_SLICE_SPACING] != NULL ){ /* get reported slice spacing */
ddd = strstr(epos[E_SLICE_SPACING],"//") ;
if( ddd != NULL ) {
if(*(ddd+2)=='\n'){ /* catch carriage returns - Jeff Gunter via DRG 3/14/2007 */
sp = 0.0; /* probably should write this as function to check on all DICOM fields*/
}
else {
sscanf( ddd+2 , "%f" , &sp ) ;
}
}
}
if( epos[E_SLICE_THICKNESS] != NULL ){ /* get reported slice thickness */
ddd = strstr(epos[E_SLICE_THICKNESS],"//") ;
if( ddd != NULL ) {
if(*(ddd+2)=='\n'){
th = 0.0;
}
else {
sscanf( ddd+2 , "%f" , &th ) ;
}
}
}
th = fabs(th) ; sp = fabs(sp) ; /* we don't use the sign */
if( stupid_ge_fix ){ /* always be stupid */
dz = sp+th ;
} else {
if( no_stupidity && sp > 0.0 ) /* 13 Jan 2004: if 'NO', then */
dz = sp ; /* always use spacing if present */
else
dz = (sp > th) ? sp : th ; /* the correct-ish DICOM way */
#define GFAC 0.99
if( !no_stupidity ){ /* unless stupidity is turned off */
if( sp > 0.0 && sp < GFAC*th ) dz = sp+th ; /* the stupid GE way again */
if( sp > 0.0 && sp < GFAC*th && nwarn < NWMAX ){
fprintf(stderr,
"++ DICOM WARNING: Slice_Spacing=%f smaller than Slice_Thickness=%f\n",
sp , th ) ;
if( nwarn == 0 )
fprintf(stderr,
"\n"
"++ Setting environment variable AFNI_SLICE_SPACING_IS_GAP ++\n"
"++ to YES will make the center-to-center slice distance ++\n"
"++ be set to Slice_Spacing+Slice_Thickness=%6.3f. ++\n"
"++ This is against the DICOM standard [attribute (0018,0088) ++\n"
"++ is defined as the center-to-center spacing between slices, ++\n"
"++ NOT as the edge-to-edge gap between slices], but it seems ++\n"
"++ to be necessary for some GE scanners. ++\n"
"++ ++\n"
"++ This correction has been made on this data: dz=%6.3f. ++\n"
"++ ++\n"
"++ Setting AFNI_SLICE_SPACING_IS_GAP to NO means that the ++\n"
"++ DICOM Slice_Spacing variable will be used for dz, replacing ++\n"
"++ the Slice_Thickness variable. This usage may be required ++\n"
"++ for some pulse sequences on Phillips scanners. ++\n"
"\n\a" ,
sp+th , dz ) ;
}
if( sp > 0.0 && sp < th && nwarn == NWMAX )
fprintf(stderr,"++ DICOM WARNING: no more Slice_Spacing messages will be printed\n") ;
nwarn++ ;
}
}
if( dz == 0.0 ) dz = 1.0 ; /* nominal dz */
return(dz);
} /*-- end of dz code, with all its stupidities --*/
/*------------------------------------------------------------------------------*/
/*! Count images in a DICOM file, if possible.
--------------------------------------------------------------------------------*/
int mri_imcount_dicom( char *fname )
{
char *ppp , *ddd ;
off_t poff ;
int plen ;
char *epos[NUM_ELIST] ;
int ii , ee , bpp , datum ;
int nx,ny,nz ;
int mosaic=0 , mos_nx,mos_ny , mos_ix,mos_iy,mos_nz ; /* 28 Oct 2002 */
Siemens_extra_info sexinfo ; /* 02 Dec 2002 */
char *sexi_start, *sexi_start2; /* KRH 25 Jul 2003 */
char *sexi_end;
ENTRY("mri_imcount_dicom") ;
if( str_sexinfo != NULL ){ free(str_sexinfo); str_sexinfo=NULL; }
if( !mri_possibly_dicom(fname) ) RETURN(0) ; /* 07 May 2003 */
/* extract the header from the file (cf. mri_dicom_hdr.[ch]) */
mri_dicom_nohex(1) ; mri_dicom_noname(1) ;
ppp = mri_dicom_header( fname ) ;
if( ppp == NULL ) RETURN(0);
/* find out where the pixel array is in the file */
mri_dicom_pxlarr( &poff , &plen ) ;
if( plen <= 0 ){ free(ppp) ; RETURN(0); }
/* check if file is actually this big */
{ unsigned int psiz , fsiz ;
fsiz = THD_filesize( fname ) ;
psiz = (unsigned int)(poff) + (unsigned int)(plen) ;
if( fsiz < psiz ){ free(ppp) ; RETURN(0); }
}
/* find positions in header of elements we care about */
for( ee=0 ; ee < NUM_ELIST ; ee++ )
epos[ee] = strstr(ppp,elist[ee]) ;
/* see if the header has the elements we absolutely need */
if( epos[E_ROWS] == NULL ||
epos[E_COLUMNS] == NULL ||
epos[E_BITS_ALLOCATED] == NULL ){ free(ppp) ; RETURN(0); }
/* check if we have 1 sample per pixel (can't deal with 3 or 4 now) */
if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
if( ii != 1 ){ free(ppp) ; RETURN(0); }
}
/* check if photometric interpretation is MONOCHROME (don't like PALETTE) */
if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
if( ddd == NULL ){ free(ppp) ; RETURN(0); }
}
/* check if we have 8, 16, or 32 bits per pixel */
ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
if( ddd == NULL ){ free(ppp); RETURN(0); }
bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
switch( bpp ){
default: free(ppp) ; RETURN(0); /* bad value */
case 8: datum = MRI_byte ; break ;
case 16: datum = MRI_short; break ;
case 32: datum = MRI_int ; break ;
}
bpp /= 8 ; /* now bytes per pixel, instead of bits */
/* get nx, ny, nz */
ddd = strstr(epos[E_ROWS],"//") ;
if( ddd == NULL ){ free(ppp) ; RETURN(0); }
nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
if( nx < 2 ){ free(ppp) ; RETURN(0); }
ddd = strstr(epos[E_COLUMNS],"//") ;
if( ddd == NULL ){ free(ppp) ; RETURN(0); }
ny = 0 ; sscanf(ddd+2,"%d",&ny) ;
if( ny < 2 ){ free(ppp) ; RETURN(0); }
nz = 0 ;
if( epos[E_NUMBER_OF_FRAMES] != NULL ){
ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
}
if( nz == 0 ) nz = plen / (bpp*nx*ny) ;
/*-- 28 Oct 2002: Check if this is a Siemens mosaic. --*/
/*-- 02 Dec 2002: Don't use Acquisition Matrix anymore;
instead, use the Siemens extra info
in epos[E_SIEMENS_2]. --*/
/*-- 24 Dec 2002: Extract str_sexinfo even if not a mosaic. --*/
if( epos[E_ID_MANUFACTURER] != NULL &&
strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
epos[E_SIEMENS_2] != NULL ){
int len=0,loc=0 , aa,bb ;
sscanf(epos[E_SIEMENS_2],"%x%x%d [%d" , &aa,&bb , &len,&loc ) ;
if( len > 0 && loc > 0 ){
FILE *fp = fopen( fname , "rb" ) ;
if( fp != NULL ){
str_sexinfo = extract_bytes_from_file( fp, (off_t)loc, (size_t)len, 1 ) ;
fclose(fp) ;
}
}
}
/* if assume_dicom_mosaic is not set, then require "MOSAIC" string */
/* 13 Mar 2006 [rickr] */
if( ( assume_dicom_mosaic ||
( epos[E_ID_IMAGE_TYPE] && strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC") ) )
&& str_sexinfo != NULL ){
/* 31 Oct 2002: extract extra Siemens info from file */
/* KRH 25 Jul 2003 if start and end markers are present for
* Siemens extra info, cut string down to those boundaries */
sexi_start = strstr(str_sexinfo, "### ASCCONV BEGIN ###");
if(sexi_start != NULL) { /* search for end after start - drg,fredtam 23 Mar 2007 */
sexi_start2 = strstr(sexi_start+21, "### ASCCONV BEGIN ###");
sexi_end = strstr(sexi_start, "### ASCCONV END ###");
if (sexi_end != NULL) {
char *sexi_tmp;
int sexi_size;
if((sexi_start2!=NULL) && (sexi_start2<sexi_end)) {
sexi_start = sexi_start2;
}
sexi_size = sexi_end - sexi_start + 19 ;
sexi_tmp = AFMALL( char, sexi_size );
memcpy(sexi_tmp,sexi_start,sexi_size);
free(str_sexinfo);
str_sexinfo = sexi_tmp;
}
}
/* end KRH 25 Jul 2003 change */
sexinfo.good = 0 ; /* start by marking it as bad */
get_siemens_extra_info( str_sexinfo , &sexinfo ) ;
if( sexinfo.good ){ /* if data is good */
if(( epos[E_ID_IMAGE_TYPE] != NULL &&
strstr(epos[E_ID_IMAGE_TYPE],"MOSAIC") != NULL ) ||
sexinfo.mosaic_num > 1 ) {
nz = sexinfo.mosaic_num ;
} /* end of if actually MOSAIC */
} /* end of if sexinfo was good */
else { /* warn if sexinfo was bad */
static int nwarn=0 ;
if( nwarn < NWMAX )
fprintf(stderr,"++ DICOM WARNING: indecipherable SIEMENS MOSAIC info (%s) in file %s\n",
elist[E_SIEMENS_2] , fname ) ;
if( nwarn == NWMAX )
fprintf(stderr,"++ DICOM NOTICE: no more SIEMENS MOSAIC info messages will be printed\n");
nwarn++ ;
}
} /* end of if str_sexinfo != NULL */
free(ppp); RETURN(nz);
}
/*--------------------------------------------------------------------------------*/
/*! Read some bytes from an open file at a given offset. Return them in a
newly malloc()-ed array. If return value is NULL, something bad happened.
----------------------------------------------------------------------------------*/
static char * extract_bytes_from_file( FILE *fp, off_t start, size_t len, int strize )
{
char *ar ;
size_t nn , ii ;
if( fp == NULL || len == 0 ) return NULL ; /* bad inputs? */
ar = AFMALL(char, len+1) ; /* make space for data */
lseek( fileno(fp) , start , SEEK_SET ) ; /* set file position */
nn = fread( ar , 1 , len , fp ) ; /* read data */
if( nn == 0 ){ free(ar); return NULL; } /* bad read? */
if( strize ){ /* convert to C string? */
for( ii=0 ; ii < nn ; ii++ ) /* scan for NULs and */
if( ar[ii] == '\0' ) ar[ii] = ' ' ; /* replace with blanks */
}
return ar ;
}
/*--------------------------------------------------------------------------------*/
/*! Parse the Siemens extra stuff for mosaic information.
Ad hoc, based on sample data and no documentation.
----------------------------------------------------------------------------------*/
static void get_siemens_extra_info( char *str , Siemens_extra_info *mi )
{
char *cpt , *dpt, *ept ;
int nn , mm , snum , last_snum=-1 ;
int have_x[2] = {0,0},
have_y[2] = {0,0},
have_z[2] = {0,0};
float x,y,z , val ;
char name[1024] ;
/*-- check for good inputs --*/
if( mi == NULL ) return ;
mi->good = 0 ;
for( snum=0 ; snum < NMOMAX ; snum++ )
mi->slice_xyz[snum][0] = mi->slice_xyz[snum][1] = mi->slice_xyz[snum][2] = -9999.9 ;
if( str == NULL || *str == '\0' ) return ;
/*-- find string that starts the slice information array --*/
/* 04 Mar 2003 reworked this section to skip "fake matches" *
* of the target string present in some mosaic files in the *
* binary section --KRH */
nn = 0;
ept = str; /* use temporary pointer instead of passed pointer to Siemens */
/* must be able to read at least 3 of the 4 parameters in slice information */
while ((nn < 3) && (strlen(ept) > 20)) { /* mod drg, fredtam */
cpt = strstr( str , "sSliceArray.asSlice[" ) ; /* 20 characters minimum */
if( cpt == NULL ) return ;
/* interpret next string into
snum = slice subscript (0,1,...)
name = variable name
val = number of RHS of '=' sign
mm = # of bytes used in scanning the above */
nn = sscanf( cpt , "sSliceArray.asSlice[%d].%1022s =%f%n" ,
&snum , name , &val , &mm ) ;
ept = cpt + 20; /* skip to end of "false match", advance to next string KRH */
}
/*-- scan for coordinates, until can't find a good string to scan --*/
#if 0
/* 25 Feb 2003 Changed logic here KRH */
have_x = have_y = have_z = 0 ;
#endif
while(1){
if( nn < 3 ) break ; /* bad conversion set */
if( snum < 0 || snum >= NMOMAX ) break ; /* slice number out of range */
/* 21 Feb 2003 rework this section to allow for missing coordinates in
* some mosaic files --KRH */
#if 0
/* assign val based on name */
if( strcmp(name,"sPosition.dSag") == 0 ){ x = val; have_x = 1; }
else if( strcmp(name,"sPosition.dCor") == 0 ){ y = val; have_y = 1; }
else if( strcmp(name,"sPosition.dTra") == 0 ){ z = val; have_z = 1; }
/* if now have all 3 coordinates, save them */
if( have_x && have_y && have_z ){
mi->slice_xyz[snum][0] = x; mi->slice_xyz[snum][1] = y; mi->slice_xyz[snum][2] = z;
last_snum = snum ;
have_x = have_y = have_z = 0 ;
}
#endif
if( strcmp(name,"sPosition.dSag") == 0 ){
mi->slice_xyz[snum][0] = val;
if (snum < 2) have_x[snum] = 1;
last_snum = snum;
} else if( strcmp(name,"sPosition.dCor") == 0 ){
mi->slice_xyz[snum][1] = val;
if (snum < 2) have_y[snum] = 1;
last_snum = snum;
} else if( strcmp(name,"sPosition.dTra") == 0 ){
mi->slice_xyz[snum][2] = val;
if (snum < 2) have_z[snum] = 1;
last_snum = snum;
}
/* skip to next slice array assignment string (which may not be a coordinate) */
dpt = cpt + mm ; /* just after 'val' */
cpt = dpt ;
while( isspace(*cpt) ) cpt++ ; /* skip over whitespace */
if( cpt-dpt > 16 ) break ; /* too much space */
if( strncmp(cpt,"sSliceArray.asSlice[",20) != 0 ) break ; /* bad next line */
/* 04 Mar 2003 moved this stuff around to allow for locating "fake matches" *
* of the target text in some mosaic files' binary sections */
/* interpret next string into
snum = slice subscript (0,1,...)
name = variable name
val = number of RHS of '=' sign
mm = # of bytes used in scanning the above */
nn = sscanf( cpt , "sSliceArray.asSlice[%d].%1022s =%f%n" ,
&snum , name , &val , &mm ) ;
}
/* if got at least 1 slice info, mark data as being good */
if( last_snum >= 0 ){
mi->good = 1 ;
if (have_x[0] && have_x[1]) mi->have_data[0] = 1;
if (have_y[0] && have_y[1]) mi->have_data[1] = 1;
if (have_z[0] && have_z[1]) mi->have_data[2] = 1;
mi->mosaic_num = last_snum+1 ;
}
return ;
}
/*--------------------------------------------------------------------------*/
/*! Test if a file is possibly a DICOM file. -- RWCox - 07 May 2003
----------------------------------------------------------------------------*/
int mri_possibly_dicom( char *fname )
{
#undef BSIZ
#define BSIZ 4096
FILE *fp ;
unsigned char buf[BSIZ] , *cpt ;
int nn , ii ;
if( fname == NULL || *fname == '\0' ) return 0 ;
fp = fopen( fname , "rb" ) ; if( fp == NULL ) return 0 ;
/* read 1st buffer */
nn = fread( buf , 1 , BSIZ , fp ) ;
if( nn < 256 ){ fclose(fp) ; return 0 ; } /* too short */
/* easy: check if has 'DICM' marker at offset 128..131 */
if( buf[128]=='D' && buf[129]=='I' && buf[130]=='C' && buf[131]=='M' ){
fclose(fp) ; return 1 ;
}
/* hard: scan file for sequence: E0 7F 10 00 (image data attribute) */
while(1){
cpt = memchr( buf, 0xe0, nn ) ; /* look for E0 */
if( cpt == NULL ){ /* skip this buffer */
nn = fread( buf , 1 , BSIZ , fp ) ; /* and get another */
if( nn < 256 ){ fclose(fp) ; return 0 ; }
continue ;
}
ii = nn - (cpt-buf) ; /* num char to end of buf */
if( ii <= 4 ){ /* too close to end of buf */
memmove( buf , cpt , ii ) ;
nn = fread( buf+ii , 1 , BSIZ-ii , fp ) ; nn += ii ;
if( nn < 256 ){ fclose(fp) ; return 0 ; }
cpt = buf ; ii = nn ;
}
/* see if we got what we want */
if( *cpt==0xe0 && *(cpt+1)==0x7f && *(cpt+2)==0x10 && *(cpt+3)==0x00 ){
fclose(fp) ; return 1 ;
}
/* no? start again at next char in buf */
memmove( buf , cpt+1 , ii-1 ) ; nn = ii-1 ;
}
}
/* clear oblique information structure */
static void
Clear_obl_info(oblique_info *obl_info)
{
int i,j;
LOAD_FVEC3(obl_info->dfpos1,0.0,0.0,0.0);
LOAD_FVEC3(obl_info->dfpos2,0.0,0.0,0.0);
LOAD_FVEC3(obl_info->del,0.0,0.0,0.0);
LOAD_FVEC3(obl_info->xvec,0.0,0.0,0.0);
LOAD_FVEC3(obl_info->yvec,0.0,0.0,0.0);
obl_info->mosaic = 0;
obl_info->mos_ix = obl_info->mos_nx = obl_info->mos_ny =
obl_info->mos_nslice = 1;
obl_info->nx = obl_info->ny = 1;
obl_info_set = 0;
/* make all elements zero flagging it hasn't been computed yet */
/* lower right corner of valid MAT44 matrix is 1.0, so this is invalid */
for(i=0;i<4;i++) {
for(j=0;j<4;j++) {
obl_info->Tr_dicom[i][j] = 0.0;
}
}
/* memset(&obl_info->Tr_dicom[0][0], 0, 16*sizeof(float)); */
}
/* fill oblique information structure */
static void
Fill_obl_info(oblique_info *obl_info, char **epos, Siemens_extra_info *siem)
{
float *xyz ; int qq ;
char *ddd;
float dx, dy, th, sp, dz;
/* float xc1, xc2, xc3, yc1, yc2, yc3, xn, yn;*/
int nx, ny, mos_ix, mos_iy;
int ii;
THD_fvec3 xc, yc;
ENTRY("Fill_obl_info");
if(obl_info_set) /* if already set all parameters for first slice */
xyz = obl_info->dfpos2.xyz; /* only need to set ImagePosition for 2nd slice */
else
xyz = obl_info->dfpos1.xyz;
if(epos[E_IMAGE_POSITION] != NULL ){ /* origin position of slice */
ddd = strstr(epos[E_IMAGE_POSITION],"//") ;
if( ddd != NULL ){
qq = sscanf(ddd+2,"%f\\%f\\%f",xyz,xyz+1,xyz+2) ;
}
}
if(obl_info_set) {
obl_info_set = 2;
EXRETURN;
}
if( epos[E_PIXEL_SPACING] != NULL ){
ddd = strstr(epos[E_PIXEL_SPACING],"//") ;
if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
if( dy == 0.0 && dx > 0.0 ) dy = dx ;
}
dz = get_dz(epos);
/* set voxel sizes */
LOAD_FVEC3(obl_info->del, dx, dy, dz);
if(epos[E_IMAGE_ORIENTATION] != NULL ){
ddd = strstr(epos[E_IMAGE_ORIENTATION],"//") ;
if( ddd != NULL ){
qq = sscanf(ddd+2,"%f\\%f\\%f\\%f\\%f\\%f",
&xc.xyz[0], &xc.xyz[1], &xc.xyz[2],
&yc.xyz[0], &yc.xyz[1], &yc.xyz[2]);
/* check if both vectors OK */
if( qq == 6 && SIZE_FVEC3(xc) > 0.0 && SIZE_FVEC3(yc) > 0.0 ){
NORMALIZE_FVEC3(xc);
NORMALIZE_FVEC3(yc);
/* if the values are close to 0 or 1 make it so */
for(ii=0;ii<3;ii++) {
if(ALMOST(xc.xyz[ii],0.0))
xc.xyz[ii] = 0.0;
if(ALMOST(xc.xyz[ii],1.0))
xc.xyz[ii] = 1.0;
if(ALMOST(xc.xyz[ii],-1.0))
xc.xyz[ii] = -1.0;
if(ALMOST(yc.xyz[ii],0.0))
yc.xyz[ii] = 0.0;
if(ALMOST(yc.xyz[ii],1.0))
yc.xyz[ii] = 1.0;
if(ALMOST(yc.xyz[ii],-1.0))
yc.xyz[ii] = -1.0;
}
obl_info->xvec = xc;
obl_info->yvec = yc;
}
}
}
obl_info_set = 1;
/* handle Siemens mosaic data */
if(siem->mosaic_num>1) {
obl_info->mosaic = 1;
obl_info->mos_sliceinfo = 0;
/* need nx, ny in mosaic and in each slice and the real number of slices*/
ddd = strstr(epos[E_COLUMNS],"//") ;
sscanf(ddd+2,"%d",&nx) ;
ddd = strstr(epos[E_ROWS],"//") ;
sscanf(ddd+2,"%d",&ny) ;
/* get siemens in the same way as in the standard mri_read_dicom above */
/* compute size of mosaic layout
as 1st integer whose square is >= # of images in mosaic */
for( mos_ix=1 ; mos_ix*mos_ix < siem->mosaic_num ; mos_ix++ ) ;
mos_iy = mos_ix ; /* number of subimages in each direction */
obl_info->mos_ix = mos_ix;
obl_info->mos_nx = nx / mos_ix ; /* sub-image dimensions*/
obl_info->mos_ny = ny / mos_iy ;
obl_info->mos_nslice = siem->mosaic_num;
obl_info->nx = nx; obl_info->ny = ny;
obl_info_set = 2;
/* check for alternate slice information */
if((siem->slice_xyz[0][0] == -9999.9) ||
(siem->slice_xyz[0][1] == -9999.9) || \
(siem->slice_xyz[0][2] == -9999.9) || \
(siem->slice_xyz[obl_info->mos_nslice-1][0] == -9999.9) ||
(siem->slice_xyz[obl_info->mos_nslice-1][1] == -9999.9) || \
(siem->slice_xyz[obl_info->mos_nslice-1][2] == -9999.9)) {
WARNING_message( \
"cannot compute alternative slice-based center for Siemens mosaic data\n");
}
else {
obl_info->slice_xyz[0][0] = siem->slice_xyz[0][0];
obl_info->slice_xyz[0][1] = siem->slice_xyz[0][1];
obl_info->slice_xyz[0][2] = siem->slice_xyz[0][2];
obl_info->slice_xyz[1][0] = siem->slice_xyz[obl_info->mos_nslice-1][0];
obl_info->slice_xyz[1][1] = siem->slice_xyz[obl_info->mos_nslice-1][1];
obl_info->slice_xyz[1][2] = siem->slice_xyz[obl_info->mos_nslice-1][2];
obl_info->mos_sliceinfo = 1;
}
}
}
/* check if data is oblique by using the vectors from the ImageOrientation field */
static int CheckObliquity(float xc1, float xc2, float xc3, float yc1, float yc2, float yc3)
{
int obliqueflag = 0;
/* any values not 1 or 0 or really close mean the data is oblique */
if ((!ALMOST(fabs(xc1),1.0) && !ALMOST(xc1,0.0)) ||
(!ALMOST(fabs(xc2),1.0) && !ALMOST(xc2,0.0)) ||
(!ALMOST(fabs(xc3),1.0) && !ALMOST(xc3,0.0)) ||
(!ALMOST(fabs(yc1),1.0) && !ALMOST(yc1,0.0)) ||
(!ALMOST(fabs(yc2),1.0) && !ALMOST(yc2,0.0)) ||
(!ALMOST(fabs(yc3),1.0) && !ALMOST(yc3,0.0)) )
obliqueflag = 1;
return(obliqueflag);
}
/* mod -16 May 2007 */
/* compute Tr transformation matrix for oblique data */
/* 16 element float array */
static float *ComputeObliquity(oblique_info *obl_info)
{
/* THD_fvec3 vec1, vec2;*/
THD_fvec3 vec3, vec4, vec5, vec6, dc1, dc2, dc3, dc4 ;
THD_fvec3 offsetxvec, offsetyvec,offsetzvec, Cm, Orgin, Cx;
/* double dotp, angle, aangle;*/
float fac;
int altsliceinfo = 0;
Siemens_extra_info *siem;
int ii,jj;
double Cxx, Cxy, Cxz;
ENTRY("ComputeObliquity");
/* compute cross product of image orientation vectors*/
vec3 = CROSS_FVEC3(obl_info->xvec, obl_info->yvec);
/* compute dfpos (difference between first and second
ImagePositionPatient fields as a vector */
vec4 = SUB_FVEC3(obl_info->dfpos2, obl_info->dfpos1);
/* scale directions by voxel sizes*/
dc1 = SCALE_FVEC3(obl_info->xvec, obl_info->del.xyz[0]);
dc2 = SCALE_FVEC3(obl_info->yvec, obl_info->del.xyz[1]);
dc3 = SCALE_FVEC3(vec3, obl_info->del.xyz[2]);
/* if not Siemens mosaic this should be enough (GE for instance)*/
if(!obl_info->mosaic) {
vec5 = NORMALIZE_FVEC3(vec4);
vec6 = NORMALIZE_FVEC3(dc3);
fac = DOT_FVEC3(vec5, vec6);
if(fac==0){
WARNING_message("Bad DICOM header - assuming oblique scaling direction!");
fac = 1;
}
else {
if((fac!=1)&&(fac!=-1)) {
WARNING_message("Image Positions do not lie in same direction as \
cross product vector - %f", fac);
}
if(fac >0) fac = 1;
else fac = -1;
}
}
else fac = 1;
/* switch direction of normal vector by factor */
dc4 = SCALE_FVEC3(dc3, fac);
#ifdef DEBUG_ON
DUMP_FVEC3("xvec", obl_info->xvec);
DUMP_FVEC3("yvec", obl_info->yvec);
DUMP_FVEC3("vec3", vec3);
DUMP_FVEC3("dfpos1", obl_info->dfpos1);
DUMP_FVEC3("dfpos2", obl_info->dfpos2);
DUMP_FVEC3("vec4", vec4);
DUMP_FVEC3("del",obl_info->del);
DUMP_FVEC3("dc1", dc1);
DUMP_FVEC3("dc2", dc2);
DUMP_FVEC3("dc3", dc3);
DUMP_FVEC3("dc4", dc4);
#endif
/* Tr = malloc(16 * sizeof(float));*/
/* *Tr = dc1.xyz[0]; *(Tr+4) = dc1.xyz[1]; *(Tr+8) = dc1.xyz[2];*/
obl_info->Tr_dicom[0][0] = dc1.xyz[0];
obl_info->Tr_dicom[1][0] = dc1.xyz[1];
obl_info->Tr_dicom[2][0] = dc1.xyz[2];
/* *(Tr+1) = dc2.xyz[0]; *(Tr+5) = dc2.xyz[1]; *(Tr+9) = dc2.xyz[2];*/
obl_info->Tr_dicom[0][1] = dc2.xyz[0];
obl_info->Tr_dicom[1][1] = dc2.xyz[1];
obl_info->Tr_dicom[2][1] = dc2.xyz[2];
/* *(Tr+2) = dc4.xyz[0]; *(Tr+6) = dc4.xyz[1]; *(Tr+10) = dc4.xyz[2];*/
obl_info->Tr_dicom[0][2] = dc4.xyz[0];
obl_info->Tr_dicom[1][2] = dc4.xyz[1];
obl_info->Tr_dicom[2][2] = dc4.xyz[2];
/* *(Tr+3) = obl_info->dfpos1.xyz[0]; *(Tr+7) = obl_info->dfpos1.xyz[1];
*(Tr+11) = obl_info->dfpos1.xyz[2];*/
obl_info->Tr_dicom[0][3] = obl_info->dfpos1.xyz[0];
obl_info->Tr_dicom[1][3] = obl_info->dfpos1.xyz[1];
obl_info->Tr_dicom[2][3] = obl_info->dfpos1.xyz[2];
/* *(Tr+12) = *(Tr+13) = *(Tr+14) = 0.0; *(Tr+15) = 1.0;*/
obl_info->Tr_dicom[3][0] = 0; obl_info->Tr_dicom[3][1] = 0;
obl_info->Tr_dicom[3][2] = 0; obl_info->Tr_dicom[3][3] = 1.0;
if(!obl_info->mosaic) {
RETURN(&(obl_info->Tr_dicom[0][0]));
}
/* for Siemens mosaic data, seen two cases */
#ifdef DEBUG_ON
printf("mos_ix %d mos_nx %d, mos_ny %d, mos_nslice %d\n" \
"nx %d ny %d\n", \
obl_info->mos_ix, obl_info->mos_nx, obl_info->mos_ny, obl_info->mos_nslice, \
obl_info->nx, obl_info->ny);
#endif
/* Siemens mosaic with full slice information */
/* will rely on ImagePosition method for now*/
/* Siemens mosaic with limited slice information - rely on ImagePosition*/
/* compute central mosaic point - in funny way */
offsetxvec = SCALE_FVEC3(dc1, ((obl_info->nx)/2.0));
offsetyvec = SCALE_FVEC3(dc2, ((obl_info->ny)/2.0));
offsetzvec = SCALE_FVEC3(dc3, ((obl_info->mos_nslice - 1.0)*fac/2.0));
Cm = ADD_FVEC3(obl_info->dfpos1, offsetxvec);
Cm = ADD_FVEC3(Cm, offsetyvec);
Cm = ADD_FVEC3(Cm, offsetzvec);
/* find origin using single slice info */
offsetxvec = SCALE_FVEC3(dc1, ((obl_info->mos_nx - 1.0)/2.0));
offsetyvec = SCALE_FVEC3(dc2, ((obl_info->mos_ny - 1.0)/2.0));
offsetzvec = SCALE_FVEC3(dc3, ((obl_info->mos_nslice - 1.0)*fac/2.0));
Orgin = SUB_FVEC3(Cm, offsetxvec);
Orgin = SUB_FVEC3(Orgin, offsetyvec);
Orgin = SUB_FVEC3(Orgin, offsetzvec);
/* check if this center of the mosaic is about the same as the slice center*/
if(obl_info->mos_sliceinfo) {
/* find center of volume as center of vector from
center of first slice to center of last slice */
Cxx = (obl_info->slice_xyz[0][0] + obl_info->slice_xyz[1][0]) / 2.0;
Cxy = (obl_info->slice_xyz[0][1] + obl_info->slice_xyz[1][1]) / 2.0;
Cxz = (obl_info->slice_xyz[0][2] + obl_info->slice_xyz[1][2]) / 2.0;
LOAD_FVEC3(Cx, Cxx, Cxy, Cxz);
/* also check if the vector from slice 0 to last slice parallel or
anti-parallel to the slice normal */
LOAD_FVEC3(obl_info->dfpos1,
(obl_info->slice_xyz[1][0] - obl_info->slice_xyz[0][0]),
(obl_info->slice_xyz[1][1] - obl_info->slice_xyz[0][1]),
(obl_info->slice_xyz[1][2] - obl_info->slice_xyz[0][2]));
vec5 = NORMALIZE_FVEC3(obl_info->dfpos1);
vec6 = NORMALIZE_FVEC3(dc3);
fac = DOT_FVEC3(vec5, vec6);
if(fac==0){
WARNING_message("Bad DICOM header - assuming oblique scaling direction!");
fac = 1;
obl_info->mos_sliceinfo = 0;
}
else {
if((fac!=1)&&(fac!=-1)) {
WARNING_message("Image Positions do not lie in same direction"
"as cross product vector - %f", fac);
}
if(fac >0) fac = 1;
else fac = -1;
obl_info->mos_sliceinfo = 1;
}
if(fac==-1) {
INFO_message("Assuming anti-parallel (left-handed coordinate system)");
}
if(!ALMOST(Cm.xyz[0], Cx.xyz[0]) ||
!ALMOST(Cm.xyz[1], Cx.xyz[1]) ||
!ALMOST(Cm.xyz[2], Cx.xyz[2])) {
WARNING_message("Slice-based center is different from mosaic center");
WARNING_message("Origin computation of obliquity may be incorrect");
DUMP_FVEC3("Mosaic Center", Cm);
DUMP_FVEC3("Slice-based Center", Cx);
}
else
INFO_message("Slice based center matches mosaic center\n");
}
DUMP_FVEC3("Mosaic Center", Cm);
DUMP_FVEC3("Origin Coordinates", Orgin);
/* update 4th column of transformation matrix with computed origin */
obl_info->Tr_dicom[0][3] = Orgin.xyz[0];
obl_info->Tr_dicom[1][3] = Orgin.xyz[1];
obl_info->Tr_dicom[2][3] = Orgin.xyz[2];
/* adjust for rounding errors by setting values close to 0 or 1 to 0 or 1 */
for(ii=0;ii<4;ii++) {
for(jj=0;jj<4;jj++) {
if(ALMOST(obl_info->Tr_dicom[ii][jj], 0.0))
obl_info->Tr_dicom[ii][jj] = 0.0;
if(ALMOST(obl_info->Tr_dicom[ii][jj], 1.0))
obl_info->Tr_dicom[ii][jj] = 1.0;
if(ALMOST(obl_info->Tr_dicom[ii][jj], -1.0))
obl_info->Tr_dicom[ii][jj] = -1.0;
}
}
#ifdef DEBUG_ON
DUMP_FVEC3("Image Position", obl_info->dfpos1);
DUMP_FVEC3("dc1", dc1);
DUMP_FVEC3("dc2", dc2);
DUMP_FVEC3("dc3", dc3);
DUMP_FVEC3("Center of Mosaic", Cm);
DUMP_FVEC3("Origin", Orgin);
#endif
RETURN(&(obl_info->Tr_dicom[0][0]));
#if 0
/* compute dot product with 0,1,0 axis - put 1 in maximum index
(closest major axis) */
vec4.xyz[0] = 0.0; vec4.xyz[1] = 0.0; vec4.xyz[2] = 0.0;
vec4.xyz[MAXINDEX_FVEC3(vec2)] = 1.0;
DUMP_FVEC3("vec4", vec4);
dotp = DOT_FVEC3(vec3, vec4);
/* compute angle as inverse cosine in degrees */
angle = -(90.0 - 180.0*acos(dotp) / PI);
aangle = abs(angle);
if(aangle<0.001)
angle = 0.0;
RETURN(angle);
#endif
}
/* externally available function to reset oblique info */
void mri_read_dicom_reset_obliquity()
{
Clear_obl_info(&obl_info);
}
/* externally available function to compute oblique transformation */
/* and store resulting matrix in 16 element float array */
void mri_read_dicom_get_obliquity(float *Tr)
{
float *fptr;
int i,j;
fptr = Tr;
if(obl_info_set) /* if oblique info is filled in, even partially */
ComputeObliquity(&obl_info); /* compute proper transformation or warn */
/* otherwise just ignore it-no warnings (for non-DICOM data) */
for(i=0;i<4;i++)
for(j=0;j<4;j++)
*fptr++ = obl_info.Tr_dicom[i][j];
/* memcpy(Tr, fptr, 16*sizeof(float));*/
return;
}
syntax highlighted by Code2HTML, v. 0.9.1