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