00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <ctype.h>
00006 #include <sys/types.h>
00007
00008 #include "mri_image.h"
00009 #include "mri_dicom_hdr.h"
00010 #include "Amalloc.h"
00011 #include "dbtrace.h"
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 extern char * mri_dicom_header( char * );
00045 extern void mri_dicom_pxlarr( off_t *, unsigned int * ) ;
00046 extern void mri_dicom_noname( int ) ;
00047 extern void mri_dicom_nohex ( int ) ;
00048
00049 static int use_DI_MRL_xcos = 0;
00050 static int use_DI_MRL_ycos = 0;
00051 static int use_DI_MRL_zcos = 0;
00052
00053 static float DI_MRL_xcos[3] = { 1.0, 0.0, 0.0 };
00054 static float DI_MRL_ycos[3] = { 0.0, 1.0, 0.0 };
00055 static float DI_MRL_zcos[3] = { 0.0, 0.0, 1.0 };
00056
00057 static float DI_MRL_xoff = 0.0;
00058 static float DI_MRL_yoff = 0.0;
00059 static float DI_MRL_zoff = 0.0;
00060 static int use_DI_MRL_xoff = 0;
00061 static int use_DI_MRL_yoff = 0;
00062 static int use_DI_MRL_zoff = 0;
00063
00064
00065
00066 char DI_MRL_orients[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
00067 float DI_MRL_tr = 0.0;
00068
00069 extern unsigned long l_THD_filesize( char * pathname );
00070
00071 void swap_twobytes( int nvals, void * data )
00072 {
00073 char * dp0 = (char *)data;
00074 char * dp1 = dp0 + 1;
00075 int c;
00076
00077 for( c = 0; c < nvals; c++, dp0 += 2, dp1 += 2 )
00078 {
00079 *dp0 ^= *dp1;
00080 *dp1 ^= *dp0;
00081 *dp0 ^= *dp1;
00082 }
00083 }
00084
00085 void swap_fourbytes( int nvals, void * data )
00086 {
00087 char * dp0 = (char *)data;
00088 char * dp1 = dp0 + 1;
00089 char * dp2 = dp0 + 2;
00090 char * dp3 = dp0 + 3;
00091 int c;
00092
00093 for( c = 0; c < nvals; c++, dp0 += 4, dp1 += 4, dp2 += 4, dp3 += 4 )
00094 {
00095 *dp0 ^= *dp3; *dp3 ^= *dp0; *dp0 ^= *dp3;
00096 *dp1 ^= *dp2; *dp2 ^= *dp1; *dp1 ^= *dp2;
00097 }
00098 }
00099
00100
00101
00102 struct dimon_stuff_t {
00103 int study, series, image;
00104 } gr_dimon_stuff;
00105
00106
00107
00108 static int LITTLE_ENDIAN_ARCHITECTURE = -1 ;
00109
00110 static void set_endianosity(void)
00111 {
00112 long val = 1;
00113
00114 LITTLE_ENDIAN_ARCHITECTURE = (*(char *)&val == 1);
00115 }
00116
00117
00118
00119 static char *elist[] = {
00120
00121 "0018 0050" ,
00122 "0018 0080" ,
00123 "0018 0088" ,
00124 "0018 1149" ,
00125
00126 "0020 0020" ,
00127 "0020 0032" ,
00128 "0020 0037" ,
00129 "0020 1041" ,
00130
00131 "0028 0002" ,
00132 "0028 0008" ,
00133 "0028 0010" ,
00134 "0028 0011" ,
00135 "0028 0030" ,
00136 "0028 0100" ,
00137 "0028 0101" ,
00138 "0028 1052" ,
00139 "0028 1053" ,
00140 "0028 1054" ,
00141 "0028 0004" ,
00142 "0028 0103" ,
00143 "0028 0102" ,
00144 "0028 1050" ,
00145 "0028 1051" ,
00146
00147 "0008 0008" ,
00148 "0008 0070" ,
00149 "0018 1310" ,
00150
00151 "0029 1010" ,
00152 "0029 1020" ,
00153
00154 "0020 0010" ,
00155 "0020 0011" ,
00156 "0020 0013" ,
00157
00158 NULL } ;
00159
00160 #define NUM_ELIST (sizeof(elist)/sizeof(char *)-1)
00161
00162 #define E_SLICE_THICKNESS 0
00163 #define E_REPETITION_TIME 1
00164 #define E_SLICE_SPACING 2
00165 #define E_FIELD_OF_VIEW 3
00166
00167 #define E_PATIENT_ORIENTATION 4
00168 #define E_IMAGE_POSITION 5
00169 #define E_IMAGE_ORIENTATION 6
00170 #define E_SLICE_LOCATION 7
00171
00172 #define E_SAMPLES_PER_PIXEL 8
00173 #define E_NUMBER_OF_FRAMES 9
00174 #define E_ROWS 10
00175 #define E_COLUMNS 11
00176 #define E_PIXEL_SPACING 12
00177 #define E_BITS_ALLOCATED 13
00178 #define E_BITS_STORED 14
00179 #define E_RESCALE_INTERCEPT 15
00180 #define E_RESCALE_SLOPE 16
00181 #define E_RESCALE_TYPE 17
00182 #define E_PHOTOMETRIC_INTERPRETATION 18
00183 #define E_PIXEL_REPRESENTATION 19
00184 #define E_HIGH_BIT 20
00185 #define E_WINDOW_CENTER 21
00186 #define E_WINDOW_WIDTH 22
00187
00188 #define E_ID_IMAGE_TYPE 23
00189 #define E_ID_MANUFACTURER 24
00190 #define E_ACQ_MATRIX 25
00191
00192 #define E_SIEMENS_1 26
00193 #define E_SIEMENS_2 27
00194
00195 #define E_RS_STUDY_NUM 28
00196 #define E_RS_SERIES_NUM 29
00197 #define E_RS_IMAGE_NUM 30
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 int mri_possibly_dicom( char *fname )
00208 {
00209 #undef BSIZ
00210 #define BSIZ 4096
00211 FILE *fp ;
00212 unsigned char buf[BSIZ] , *cpt ;
00213 int nn , ii ;
00214
00215 if( fname == NULL || *fname == '\0' ) return 0 ;
00216 fp = fopen( fname , "rb" ) ; if( fp == NULL ) return 0 ;
00217
00218
00219
00220 nn = fread( buf , 1 , BSIZ , fp ) ;
00221 if( nn < 256 ){ fclose(fp) ; return 0 ; }
00222
00223
00224
00225 if( buf[128]=='D' && buf[129]=='I' && buf[130]=='C' && buf[131]=='M' ){
00226 fclose(fp) ; return 1 ;
00227 }
00228
00229
00230
00231 while(1){
00232
00233 cpt = memchr( buf, 0xe0, nn ) ;
00234
00235 if( cpt == NULL ){
00236 nn = fread( buf , 1 , BSIZ , fp ) ;
00237 if( nn < 256 ){ fclose(fp) ; return 0 ; }
00238 continue ;
00239 }
00240
00241 ii = nn - (cpt-buf) ;
00242 if( ii <= 4 ){
00243 memmove( buf , cpt , ii ) ;
00244 nn = fread( buf+ii , 1 , BSIZ-ii , fp ) ; nn += ii ;
00245 if( nn < 256 ){ fclose(fp) ; return 0 ; }
00246 cpt = buf ; ii = nn ;
00247 }
00248
00249
00250
00251 if( *cpt==0xe0 && *(cpt+1)==0x7f && *(cpt+2)==0x10 && *(cpt+3)==0x00 ){
00252 fclose(fp) ; return 1 ;
00253 }
00254
00255
00256
00257 memmove( buf , cpt+1 , ii-1 ) ; nn = ii-1 ;
00258 }
00259 }
00260
00261
00262
00263
00264 MRI_IMAGE * r_mri_read_dicom( char *fname, int debug, void ** data )
00265 {
00266 char *ppp , *ddd ;
00267 off_t poff ;
00268 unsigned int plen ;
00269 char *epos[NUM_ELIST] ;
00270 int ii,jj , ee , bpp , datum ;
00271 int nx,ny,nz , swap ;
00272 float dx,dy,dz,dt ;
00273 MRI_IMAGE *im ;
00274 void *iar ;
00275 FILE *fp ;
00276 int have_orients=0 ;
00277 int ior,jor,kor ;
00278 static int nzoff=0 ;
00279
00280 float dxx,dyy,dzz ;
00281
00282 char *eee ;
00283 float rescale_slope=0.0 , rescale_inter=0.0 ;
00284
00285 ENTRY("mri_read_dicom") ;
00286
00287 if( !mri_possibly_dicom(fname) )
00288 {
00289 if(debug > 2) fprintf(stderr,"** MRD: mri_possibly_dicom() failure\n");
00290 RETURN(NULL) ;
00291 }
00292
00293
00294
00295
00296
00297 mri_dicom_nohex(1) ;
00298 mri_dicom_noname(1) ;
00299 ppp = mri_dicom_header( fname ) ;
00300 if( ppp == NULL )
00301 {
00302 if(debug > 2) fprintf(stderr,"** MRD: mri_dicom_hdr() failure\n");
00303 RETURN(NULL);
00304 }
00305
00306
00307
00308 mri_dicom_pxlarr( &poff , &plen ) ;
00309 if( plen <= 0 ){
00310 if(debug > 2) fprintf(stderr,"** MRD: bad plen, %u\n", plen);
00311 free(ppp) ;
00312 RETURN(NULL);
00313 }
00314 if( debug > 3 )
00315 fprintf(stderr,"-d dicom: poff, plen = %d, %d\n",(int)poff,(int)plen);
00316
00317
00318
00319 { unsigned int psiz , fsiz ;
00320 fsiz = l_THD_filesize( fname ) ;
00321 psiz = (unsigned int)(poff) + plen ;
00322 if( fsiz < psiz ){
00323 if(debug > 2) fprintf(stderr,"** MRD: fsiz < psiz (%d,%d)\n",fsiz,psiz);
00324 free(ppp) ;
00325 RETURN(NULL);
00326 }
00327 }
00328
00329
00330
00331 for( ee=0 ; ee < NUM_ELIST ; ee++ )
00332 epos[ee] = strstr(ppp,elist[ee]) ;
00333
00334
00335
00336 if( epos[E_ROWS] == NULL ||
00337 epos[E_COLUMNS] == NULL ||
00338 epos[E_BITS_ALLOCATED] == NULL ){
00339 if(debug > 2)
00340 fprintf(stderr,"** MRD: missing ROWS, COLS or BITS (%p,%p,%p)\n",
00341 epos[E_ROWS], epos[E_COLUMNS], epos[E_BITS_ALLOCATED]);
00342 free(ppp) ;
00343 RETURN(NULL);
00344 }
00345
00346
00347
00348 if( epos[E_SAMPLES_PER_PIXEL] != NULL ){
00349 ddd = strstr(epos[E_SAMPLES_PER_PIXEL],"//") ;
00350 if( ddd == NULL ){
00351 if(debug > 2) fprintf(stderr,"** MRD: missing E_SAMPLES_PER_PIXEL\n");
00352 free(ppp) ;
00353 RETURN(NULL);
00354 }
00355 ii = 0 ; sscanf(ddd+2,"%d",&ii) ;
00356 if( ii != 1 ){
00357 if(debug > 2) fprintf(stderr,"** MRD: SAM_PER_PIX != 1, %d\n",ii);
00358 free(ppp) ;
00359 RETURN(NULL);
00360 }
00361 }
00362
00363
00364
00365 if( epos[E_PHOTOMETRIC_INTERPRETATION] != NULL ){
00366 ddd = strstr(epos[E_PHOTOMETRIC_INTERPRETATION],"MONOCHROME") ;
00367 if( ddd == NULL ){
00368 if(debug > 2) fprintf(stderr,"** MRD: photometric not monochrome\n");
00369 free(ppp) ;
00370 RETURN(NULL);
00371 }
00372 }
00373
00374
00375
00376 ddd = strstr(epos[E_BITS_ALLOCATED],"//") ;
00377 if( ddd == NULL ){
00378 if(debug > 2) fprintf(stderr,"** MRD: missing BITS_ALLOCATED\n");
00379 free(ppp);
00380 RETURN(NULL);
00381 }
00382 bpp = 0 ; sscanf(ddd+2,"%d",&bpp) ;
00383 switch( bpp ){
00384 default:
00385 if(debug > 2) fprintf(stderr,"** MRD: bad bpp value, %d\n",bpp);
00386 free(ppp); RETURN(NULL);
00387 case 8: datum = MRI_byte ; break ;
00388 case 16: datum = MRI_short; break ;
00389 case 32: datum = MRI_int ; break ;
00390 }
00391 bpp /= 8 ;
00392
00393 if( debug > 3 ) fprintf(stderr,"-d dicom: datum %d\n",datum);
00394
00395
00396
00397
00398
00399 #define NWMAX 2
00400 if( epos[E_BITS_STORED] != NULL && epos[E_HIGH_BIT] != NULL ){
00401 int bs=0 , hb=0 ;
00402 ddd = strstr(epos[E_BITS_STORED],"//") ; sscanf(ddd+2,"%d",&bs) ;
00403 ddd = strstr(epos[E_HIGH_BIT],"//") ; sscanf(ddd+2,"%d",&hb) ;
00404 if( bs != hb+1 ){
00405 static int nwarn=0 ;
00406 if( nwarn < NWMAX )
00407 fprintf(stderr,
00408 "++ DICOM WARNING: file %s has Bits_Stored=%d and High_Bit=%d\n",
00409 fname,bs,hb) ;
00410 if( nwarn == NWMAX )
00411 fprintf(stderr,"++ DICOM WARNING: no more Bits_Stored messages will be printed\n") ;
00412 nwarn++ ;
00413 }
00414 }
00415
00416
00417
00418
00419 if( epos[E_RESCALE_INTERCEPT] != NULL && epos[E_RESCALE_SLOPE] != NULL ){
00420 if( debug > 0 ){
00421 static int nwarn=0 ;
00422 if( nwarn == 0 ){
00423 fprintf(stderr,"++ DICOM file has Rescale tags, enforcing them...\n");
00424 fprintf(stderr," (no more Rescale tags messages will be printed)\n");
00425 }
00426 nwarn++ ;
00427 }
00428 ddd = strstr(epos[E_RESCALE_INTERCEPT],"//") ;
00429 sscanf(ddd+2,"%f",&rescale_inter) ;
00430 ddd = strstr(epos[E_RESCALE_SLOPE ],"//") ;
00431 sscanf(ddd+2,"%f",&rescale_slope) ;
00432 }
00433
00434
00435
00436
00437
00438 ddd = strstr(epos[E_ROWS],"//") ;
00439 if( ddd == NULL ){
00440 if(debug > 2) fprintf(stderr,"** MRD: missing E_ROWS\n");
00441 free(ppp) ;
00442 RETURN(NULL);
00443 }
00444 ny = 0 ; sscanf(ddd+2,"%d",&ny) ;
00445 if( ny < 2 ){
00446 if(debug > 2) fprintf(stderr,"** MRD: bad ny = %d\n",ny);
00447 free(ppp) ;
00448 RETURN(NULL);
00449 }
00450
00451 ddd = strstr(epos[E_COLUMNS],"//") ;
00452 if( ddd == NULL ){
00453 if(debug > 2) fprintf(stderr,"** MRD: missing E_COLUMNS\n");
00454 free(ppp) ;
00455 RETURN(NULL);
00456 }
00457 nx = 0 ; sscanf(ddd+2,"%d",&nx) ;
00458 if( nx < 2 ){
00459 if(debug > 2) fprintf(stderr,"** MRD: bad nx = %d\n",nx);
00460 free(ppp) ;
00461 RETURN(NULL);
00462 }
00463
00464
00465
00466 nz = 0 ;
00467 if( epos[E_NUMBER_OF_FRAMES] != NULL ){
00468 ddd = strstr(epos[E_NUMBER_OF_FRAMES],"//") ;
00469 if( ddd != NULL ) sscanf(ddd+2,"%d",&nz) ;
00470 if(debug>3) fprintf(stderr,"-d number of frames = %d\n",nz);
00471 }
00472
00473
00474
00475 if( nz == 0 ) nz = plen / (bpp*nx*ny) ;
00476 if( nz == 0 ){
00477 if(debug > 2) fprintf(stderr,"** MRD: bad nz = %d\n", nz);
00478 free(ppp) ;
00479 RETURN(NULL);
00480 }
00481
00482 if( nz != 1 ){
00483 fprintf(stderr,"** MRD: nz = %d, plen,bpp,nx,ny = %d,%d,%d,%d\n",
00484 nz, plen, bpp, nx, ny);
00485 free(ppp);
00486 RETURN(NULL);
00487 }
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 if( epos[E_ID_MANUFACTURER] != NULL &&
00498 strstr(epos[E_ID_MANUFACTURER],"SIEMENS") != NULL &&
00499 epos[E_SIEMENS_2] != NULL ){
00500 if( debug > 3 ) fprintf(stderr,"-d ignoring sexinfo\n");
00501 }
00502
00503
00504
00505 dx = dy = dz = dt = 0.0 ;
00506
00507
00508
00509 if( epos[E_PIXEL_SPACING] != NULL ){
00510 ddd = strstr(epos[E_PIXEL_SPACING],"//") ;
00511 if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00512 if( debug > 3 )
00513 fprintf(stderr,"-d dicom: read dx,dy from PIXEL_SP: %f, %f\n", dx, dy);
00514 if( dy == 0.0 && dx > 0.0 ) dy = dx ;
00515 }
00516 if( dx == 0.0 && epos[E_FIELD_OF_VIEW] != NULL ){
00517 ddd = strstr(epos[E_FIELD_OF_VIEW],"//") ;
00518 if( ddd != NULL ) sscanf( ddd+2 , "%f\\%f" , &dx , &dy ) ;
00519 if( debug > 3 )
00520 fprintf(stderr,"-d dicom: read dx,dy from FOV: %f, %f\n", dx, dy);
00521 if( dx > 0.0 ){
00522 if( dy == 0.0 ) dy = dx ;
00523 dx /= nx ; dy /= ny ;
00524 }
00525 }
00526
00527
00528
00529
00530 { int stupid_ge_fix , no_stupidity ;
00531 float sp=0.0 , th=0.0 ;
00532 static int nwarn=0 ;
00533
00534 eee = getenv("AFNI_SLICE_SPACING_IS_GAP") ;
00535 stupid_ge_fix = (eee != NULL && (*eee=='Y' || *eee=='y') );
00536 no_stupidity = (eee != NULL && (*eee=='N' || *eee=='n') );
00537
00538 if( epos[E_SLICE_SPACING] != NULL ){
00539 ddd = strstr(epos[E_SLICE_SPACING],"//") ;
00540 if( ddd != NULL ) sscanf( ddd+2 , "%f" , &sp ) ;
00541 }
00542 if( epos[E_SLICE_THICKNESS] != NULL ){
00543 ddd = strstr(epos[E_SLICE_THICKNESS],"//") ;
00544 if( ddd != NULL ) sscanf( ddd+2 , "%f" , &th ) ;
00545 }
00546
00547 if(debug>3) fprintf(stderr,"-d dicom: thick, spacing = %f, %f\n", th, sp);
00548
00549 th = fabs(th) ; sp = fabs(sp) ;
00550
00551 if( stupid_ge_fix ){
00552 dz = sp+th ;
00553 } else {
00554
00555 if( no_stupidity && sp > 0.0 )
00556 dz = sp ;
00557 else
00558 dz = (sp > th) ? sp : th ;
00559
00560 #define GFAC 0.99
00561
00562 if( !no_stupidity ){
00563 if( sp > 0.0 && sp < GFAC*th ) dz = sp+th;
00564
00565 if( sp > 0.0 && sp < GFAC*th && nwarn < NWMAX ){
00566 fprintf(stderr,
00567 "++ DICOM WARNING: file %s has Slice_Spacing=%f smaller than Slice_Thickness=%f\n",
00568 fname , sp , th ) ;
00569 if( nwarn == 0 )
00570 fprintf(stderr,
00571 "\n"
00572 "++ Setting environment variable AFNI_SLICE_SPACING_IS_GAP ++\n"
00573 "++ to YES will make the center-to-center slice distance ++\n"
00574 "++ be set to Slice_Spacing+Slice_Thickness=%6.3f. ++\n"
00575 "++ This is against the DICOM standard [attribute (0018,0088) ++\n"
00576 "++ is defined as the center-to-center spacing between slices, ++\n"
00577 "++ NOT as the edge-to-edge gap between slices], but it seems ++\n"
00578 "++ to be necessary for some GE scanners. ++\n"
00579 "++ ++\n"
00580 "++ This correction has been made on this data: dz=%6.3f. ++\n"
00581 "++ ++\n"
00582 "++ Setting AFNI_SLICE_SPACING_IS_GAP to NO means that the ++\n"
00583 "++ DICOM Slice_Spacing variable will be used for dz, replacing ++\n"
00584 "++ the Slice_Thickness variable. This usage may be required ++\n"
00585 "++ for some pulse sequences on Phillips scanners. ++\n"
00586 "\n\a" ,
00587 sp+th , dz ) ;
00588 }
00589 if( sp > 0.0 && sp < th && nwarn == NWMAX )
00590 fprintf(stderr,"++ DICOM WARNING: no more Slice_Spacing messages will be printed\n") ;
00591 nwarn++ ;
00592 }
00593 }
00594 if( dz == 0.0 && dx != 0.0 ) dz = 1.0 ;
00595
00596 }
00597
00598 if(debug > 3) fprintf(stderr,"-d dicom: using dxyz = %f, %f, %f\n",dx,dy,dz);
00599
00600
00601
00602 if( epos[E_REPETITION_TIME] != NULL ){
00603 ddd = strstr(epos[E_REPETITION_TIME],"//") ;
00604 if( ddd != NULL ){
00605 sscanf( ddd+2 , "%f" , &dt ) ;
00606 dt *= 0.001 ;
00607 if(debug > 3) fprintf(stderr,"-d dicom: rep. time dt = %f sec.\n",dt);
00608 }
00609 }
00610
00611
00612
00613 if( epos[E_RS_STUDY_NUM] != NULL ){
00614 ddd = strstr(epos[E_RS_STUDY_NUM],"//") ;
00615 if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.study);
00616 }
00617
00618 if( epos[E_RS_SERIES_NUM] != NULL ){
00619 ddd = strstr(epos[E_RS_SERIES_NUM],"//") ;
00620 if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.series);
00621 }
00622
00623 if( epos[E_RS_IMAGE_NUM] != NULL ){
00624 ddd = strstr(epos[E_RS_IMAGE_NUM],"//") ;
00625 if( ddd != NULL ) sscanf( ddd+2 , "%d" , &gr_dimon_stuff.image);
00626 }
00627
00628
00629
00630 fp = fopen( fname , "rb" ) ;
00631 if( fp == NULL ){
00632 if(debug > 2) fprintf(stderr,"** MRD: failed to open file '%s'\n",fname);
00633 free(ppp) ;
00634 RETURN(NULL);
00635 }
00636 lseek( fileno(fp) , poff , SEEK_SET ) ;
00637
00638
00639
00640 set_endianosity() ; swap = !LITTLE_ENDIAN_ARCHITECTURE ;
00641
00642
00643
00644 {
00645 im = (MRI_IMAGE *)calloc(1, sizeof(MRI_IMAGE));
00646 if( !im ){
00647 fprintf(stderr,"** MRD: im malloc failure!\n");
00648 free(ppp);
00649 RETURN(NULL);
00650 }
00651 im->nx = nx; im->ny = ny;
00652 im->nxy = nx * ny;
00653 im->nz = im->nt = im->nu = im->nv = im->nw = 1;
00654 im->nxyz = im->nxyzt = im->nvox = im->nxy;
00655 im->kind = datum;
00656 im->dx = im->dy = im->dz = im->dt = im->du = im->dv = 1.0;
00657 im->dw = -666.0;
00658 im->pixel_size = bpp;
00659 if( data ){
00660 if( ! *data ){
00661 *data = calloc(im->nvox, im->pixel_size);
00662 if( debug > 3 )
00663 fprintf(stderr,"+d dicom: image data alloc: %d x %d bytes\n",
00664 im->nvox, im->pixel_size);
00665 }
00666 if( ! *data ){
00667 fprintf(stderr,"** MRD: image data alloc failure\n");
00668 free(ppp);
00669 RETURN(NULL);
00670 }
00671 }
00672 }
00673
00674 if( dx > 0.0 && dy > 0.0 && dz > 0.0 ){
00675 im->dx = dx; im->dy = dy; im->dz = dz; im->dw = 1.0;
00676 }
00677 if( dt > 0.0 ) im->dt = dt ;
00678
00679 if( !data ) fclose(fp) ;
00680 else{
00681 iar = *data;
00682 fread( iar , bpp , nx*ny , fp ) ;
00683
00684 if( swap ){
00685 switch( im->pixel_size ){
00686 default: break ;
00687 case 2: swap_twobytes ( im->nvox, iar ) ; break ;
00688 case 4: swap_fourbytes( im->nvox, iar ) ; break ;
00689 case 8: swap_fourbytes( 2*im->nvox, iar ) ; break ;
00690 }
00691 im->was_swapped = 1 ;
00692 }
00693
00694
00695
00696 fclose(fp) ;
00697
00698
00699
00700 if( rescale_slope > 0.0 ){
00701 for( ii=0 ; ii < 1 ; ii++ ){
00702 switch( im->kind ){
00703 case MRI_byte:{
00704 byte *ar = (byte *)*data;
00705 for( jj=0 ; jj < im->nvox ; jj++ )
00706 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00707 }
00708 break ;
00709
00710 case MRI_short:{
00711 short *ar = (short *)*data;
00712 for( jj=0 ; jj < im->nvox ; jj++ )
00713 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00714 }
00715 break ;
00716
00717 case MRI_int:{
00718 int *ar = (int *)*data;
00719 for( jj=0 ; jj < im->nvox ; jj++ )
00720 ar[jj] = rescale_slope*ar[jj] + rescale_inter ;
00721 }
00722 break ;
00723 default:{
00724 fprintf(stderr,"** MRD: bad kind case (%d) for rescale_slope\n",
00725 im->kind);
00726 free(ppp); free(*data); *data = NULL;
00727 RETURN(NULL);
00728 }
00729 break ;
00730 }
00731 }
00732 }
00733 }
00734
00735
00736
00737
00738
00739 if( dt > 0.0 && DI_MRL_tr <= 0.0 ) DI_MRL_tr = dt ;
00740
00741
00742
00743 if( epos[E_IMAGE_ORIENTATION] != NULL ){
00744
00745
00746 ddd = strstr(epos[E_IMAGE_ORIENTATION],"//") ;
00747 if( ddd != NULL ){
00748 float xc1=0.0,xc2=0.0,xc3=0.0 , yc1=0.0,yc2=0.0,yc3=0.0 ;
00749 float xn,yn ; int qq ;
00750 qq=sscanf(ddd+2,"%f\\%f\\%f\\%f\\%f\\%f",&xc1,&xc2,&xc3,&yc1,&yc2,&yc3);
00751 xn = sqrt( xc1*xc1 + xc2*xc2 + xc3*xc3 ) ;
00752 yn = sqrt( yc1*yc1 + yc2*yc2 + yc3*yc3 ) ;
00753 if( qq == 6 && xn > 0.0 && yn > 0.0 ){
00754
00755 xc1 /= xn ; xc2 /= xn ; xc3 /= xn ;
00756 yc1 /= yn ; yc2 /= yn ; yc3 /= yn ;
00757
00758 if( !use_DI_MRL_xcos ){
00759 DI_MRL_xcos[0] = xc1 ; DI_MRL_xcos[1] = xc2 ;
00760 DI_MRL_xcos[2] = xc3 ; use_DI_MRL_xcos = 1 ;
00761 }
00762
00763 if( !use_DI_MRL_ycos ){
00764 DI_MRL_ycos[0] = yc1 ; DI_MRL_ycos[1] = yc2 ;
00765 DI_MRL_ycos[2] = yc3 ; use_DI_MRL_ycos = 1 ;
00766 }
00767
00768
00769
00770
00771
00772 dxx = fabs(xc1) ; ior = 1 ;
00773 dyy = fabs(xc2) ; if( dyy > dxx ){ ior=2; dxx=dyy; }
00774 dzz = fabs(xc3) ; if( dzz > dxx ){ ior=3; }
00775 dxx = DI_MRL_xcos[ior-1] ; if( dxx < 0. ) ior = -ior;
00776
00777 if( DI_MRL_orients[0] == '\0' ){
00778 switch( ior ){
00779 case -1: DI_MRL_orients[0] = 'L'; DI_MRL_orients[1] = 'R'; break;
00780 case 1: DI_MRL_orients[0] = 'R'; DI_MRL_orients[1] = 'L'; break;
00781 case -2: DI_MRL_orients[0] = 'P'; DI_MRL_orients[1] = 'A'; break;
00782 case 2: DI_MRL_orients[0] = 'A'; DI_MRL_orients[1] = 'P'; break;
00783 case 3: DI_MRL_orients[0] = 'I'; DI_MRL_orients[1] = 'S'; break;
00784 case -3: DI_MRL_orients[0] = 'S'; DI_MRL_orients[1] = 'I'; break;
00785 default: DI_MRL_orients[0] ='\0'; DI_MRL_orients[1] ='\0'; break;
00786 }
00787 }
00788
00789
00790
00791
00792
00793 dxx = fabs(yc1) ; jor = 1 ;
00794 dyy = fabs(yc2) ; if( dyy > dxx ){ jor=2; dxx=dyy; }
00795 dzz = fabs(yc3) ; if( dzz > dxx ){ jor=3; }
00796 dyy = DI_MRL_ycos[jor-1] ; if( dyy < 0. ) jor = -jor;
00797 if( DI_MRL_orients[2] == '\0' ){
00798 switch( jor ){
00799 case -1: DI_MRL_orients[2] = 'L'; DI_MRL_orients[3] = 'R'; break;
00800 case 1: DI_MRL_orients[2] = 'R'; DI_MRL_orients[3] = 'L'; break;
00801 case -2: DI_MRL_orients[2] = 'P'; DI_MRL_orients[3] = 'A'; break;
00802 case 2: DI_MRL_orients[2] = 'A'; DI_MRL_orients[3] = 'P'; break;
00803 case 3: DI_MRL_orients[2] = 'I'; DI_MRL_orients[3] = 'S'; break;
00804 case -3: DI_MRL_orients[2] = 'S'; DI_MRL_orients[3] = 'I'; break;
00805 default: DI_MRL_orients[2] ='\0'; DI_MRL_orients[3] ='\0'; break;
00806 }
00807 }
00808
00809 DI_MRL_orients[6] = '\0' ;
00810
00811 kor = 6 - abs(ior)-abs(jor) ;
00812
00813 have_orients = 1 ;
00814
00815 if( debug > 3 )
00816 fprintf(stderr,"-d dicom: DI_MRL_orients = %s, [ijk]or = %d,%d,%d\n",
00817 DI_MRL_orients, ior, jor, kor);
00818 }
00819 else if( debug > 3 )
00820 fprintf(stderr,"-d dicom: bad orient vectors qq = %d, xn,yn = %f,%f\n",
00821 qq, xn, yn);
00822 }
00823
00824 } else if( epos[E_PATIENT_ORIENTATION] != NULL ){
00825
00826 ddd = strstr(epos[E_PATIENT_ORIENTATION],"//") ;
00827 if( ddd != NULL ){
00828 char xc='\0' , yc='\0' ;
00829 sscanf(ddd+2,"%c\\%c",&xc,&yc) ;
00830 switch( toupper(xc) ){
00831 case 'L': DI_MRL_orients[0]='L'; DI_MRL_orients[1]='R'; ior=-1; break;
00832 case 'R': DI_MRL_orients[0]='R'; DI_MRL_orients[1]='L'; ior= 1; break;
00833 case 'P': DI_MRL_orients[0]='P'; DI_MRL_orients[1]='A'; ior=-2; break;
00834 case 'A': DI_MRL_orients[0]='A'; DI_MRL_orients[1]='P'; ior= 2; break;
00835 case 'F': DI_MRL_orients[0]='I'; DI_MRL_orients[1]='S'; ior= 3; break;
00836
00837 case 'H': DI_MRL_orients[0]='S'; DI_MRL_orients[1]='I'; ior=-3; break;
00838
00839 default: DI_MRL_orients[0]='\0';DI_MRL_orients[1]='\0'; ior= 0; break;
00840 }
00841 switch( toupper(yc) ){
00842 case 'L': DI_MRL_orients[2]='L'; DI_MRL_orients[3]='R'; jor=-1; break;
00843 case 'R': DI_MRL_orients[2]='R'; DI_MRL_orients[3]='L'; jor= 1; break;
00844 case 'P': DI_MRL_orients[2]='P'; DI_MRL_orients[3]='A'; jor=-2; break;
00845 case 'A': DI_MRL_orients[2]='A'; DI_MRL_orients[3]='P'; jor= 2; break;
00846 case 'F': DI_MRL_orients[2]='I'; DI_MRL_orients[3]='S'; jor= 3; break;
00847 case 'H': DI_MRL_orients[2]='S'; DI_MRL_orients[3]='I'; jor=-3; break;
00848 default: DI_MRL_orients[2]='\0';DI_MRL_orients[3]='\0'; jor= 0; break;
00849 }
00850 DI_MRL_orients[6]='\0' ;
00851 kor = 6 - abs(ior)-abs(jor) ;
00852 have_orients = (ior != 0 && jor != 0) ;
00853
00854 if( debug > 3 )
00855 fprintf(stderr,"-d dicom: DI_MRL_orients P = %s, [ijk]or = %d,%d,%d\n",
00856 DI_MRL_orients, ior, jor, kor);
00857 }
00858
00859 }
00860
00861 if( !have_orients ){
00862 fprintf(stderr,"** MRD: failed to determine dicom image orientation\n");
00863 free(ppp); free(im);
00864 if( data ){ free(*data); *data = NULL; }
00865 RETURN(NULL);
00866 }
00867
00868
00869
00870
00871
00872
00873 if( epos[E_IMAGE_POSITION] == NULL ){
00874 fprintf(stderr,"** MRD: missing image position\n");
00875 free(ppp); free(im);
00876 if( data ){ free(*data); *data = NULL; }
00877 RETURN(NULL);
00878 }
00879
00880 ddd = strstr(epos[E_IMAGE_POSITION],"//") ;
00881 if( ddd == NULL ){
00882 fprintf(stderr,"** MRD: missing IMAGE_POSITION\n");
00883 free(ppp); free(im);
00884 if( data ){ free(*data); *data = NULL; }
00885 RETURN(NULL);
00886 }
00887
00888 {
00889 float xyz[3] ; int qq ;
00890 qq = sscanf(ddd+2,"%f\\%f\\%f",xyz,xyz+1,xyz+2) ;
00891 if( qq != 3 ){
00892 fprintf(stderr,"** MRD: failed to read IMAGE_POSITION\n");
00893 free(ppp); free(im);
00894 if( data ){ free(*data); *data = NULL; }
00895 RETURN(NULL);
00896 }
00897
00898 im->xo = xyz[abs(ior)-1];
00899 im->yo = xyz[abs(jor)-1];
00900 im->zo = xyz[abs(kor)-1];
00901
00902 if( debug > 3 )
00903 fprintf(stderr,"-d dicom: read RAI image position %f, %f, %f\n"
00904 " and dset image position %f, %f, %f\n",
00905 xyz[0], xyz[1], xyz[2], im->xo, im->yo, im->zo );
00906
00907
00908 {
00909 static float zoff ;
00910 float zz = xyz[kor-1] ;
00911
00912 if( nzoff == 0 ){
00913
00914 zoff = zz ;
00915
00916
00917
00918 if( DI_MRL_orients[4] == '\0' ){
00919 switch( kor ){
00920 case 1: DI_MRL_orients[4] = 'L'; DI_MRL_orients[5] = 'R'; break;
00921 case 2: DI_MRL_orients[4] = 'P'; DI_MRL_orients[5] = 'A'; break;
00922 case 3: DI_MRL_orients[4] = 'I'; DI_MRL_orients[5] = 'S'; break;
00923 default: DI_MRL_orients[4] ='\0'; DI_MRL_orients[5] ='\0'; break;
00924 }
00925 DI_MRL_orients[6] = '\0' ;
00926 }
00927
00928
00929
00930 qq = abs(ior) ;
00931 DI_MRL_xoff = xyz[qq-1] ; use_DI_MRL_xoff = 1 ;
00932 if( ior > 0 ) DI_MRL_xoff = -DI_MRL_xoff ;
00933
00934 qq = abs(jor) ;
00935 DI_MRL_yoff = xyz[qq-1] ; use_DI_MRL_yoff = 1 ;
00936 if( jor > 0 ) DI_MRL_yoff = -DI_MRL_yoff ;
00937
00938 } else if( nzoff == 1 && !use_DI_MRL_zoff ){
00939
00940 float qoff = zz - zoff ;
00941 if( qoff < 0 ) kor = -kor ;
00942 #if 0
00943 fprintf(stderr," nzoff=1 kor=%d qoff=%f\n",kor,qoff) ;
00944 #endif
00945 switch( kor ){
00946 case 1: DI_MRL_orients[4] = 'R'; DI_MRL_orients[5] = 'L'; break;
00947 case -1: DI_MRL_orients[4] = 'L'; DI_MRL_orients[5] = 'R'; break;
00948 case 2: DI_MRL_orients[4] = 'A'; DI_MRL_orients[5] = 'P'; break;
00949 case -2: DI_MRL_orients[4] = 'P'; DI_MRL_orients[5] = 'A'; break;
00950 case 3: DI_MRL_orients[4] = 'I'; DI_MRL_orients[5] = 'S'; break;
00951 case -3: DI_MRL_orients[4] = 'S'; DI_MRL_orients[5] = 'I'; break;
00952 default: DI_MRL_orients[4] ='\0'; DI_MRL_orients[5] ='\0'; break;
00953 }
00954 DI_MRL_orients[6] = '\0' ;
00955
00956
00957
00958
00959
00960 DI_MRL_zoff = zoff ; use_DI_MRL_zoff = 1 ;
00961 if( kor > 0 ) DI_MRL_zoff = -DI_MRL_zoff ;
00962 }
00963 nzoff++ ;
00964 }
00965 }
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977 {
00978 static int nwarn = 0;
00979
00980 if( ( epos[E_SLICE_LOCATION] == NULL) ||
00981 ( (ddd = strstr(epos[E_SLICE_LOCATION],"//")) == NULL ) )
00982 {
00983 if( nwarn == 0 && debug > 1 )
00984 fprintf(stderr,"-d dimon: missing SLICE_LOCATION, continuing...\n");
00985 nwarn++;
00986 } else {
00987
00988 float zz ; int qq;
00989 qq = sscanf(ddd+2,"%f",&zz) ;
00990 if( qq != 1 ){
00991 if( !nwarn )fprintf(stderr,"** failed to extract SLICE_LOCATION\n");
00992 nwarn++;
00993 } else {
00994
00995 if( zz * im->zo < 0.0 ){
00996 if( (nwarn == 0) && (debug > 3) )
00997 fprintf(stderr,"-d image and slice loc diff in sign: %f, %f\n",
00998 im->zo, zz);
00999 nwarn++;
01000 zz = -zz;
01001 }
01002
01003 if( fabs(zz - im->zo) > 0.1 ){
01004 fprintf(stderr,
01005 "** MRD: IMAGE_LOCATION and SLICE_LOCATION disagree!\n"
01006 " z coord from IL = %f, from SL = %f\n", im->zo,zz);
01007 free(ppp); free(im);
01008 if( data ){ free(*data); *data = NULL; }
01009 RETURN(NULL);
01010 }
01011
01012 if( debug > 3 )
01013 fprintf(stderr,"-d dicom: using slice location %f (zo = %f)\n",
01014 zz, im->zo );
01015 im->zo = zz;
01016 }
01017 }
01018 }
01019
01020 free(ppp);
01021
01022 RETURN( im );
01023 }