80 : DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
81 fCompression(0), fNFiles(0), fRows(0), fColumns(0),
82 fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
83 fPixelSpacingX(0.), fPixelSpacingY(0.),
84 fSliceThickness(0.), fSliceLocation(0.),
85 fRescaleIntercept(0), fRescaleSlope(0),
86 fLittleEndian(true), fImplicitEndian(false),
87 fPixelRepresentation(0), fNbrequali(0),
88 fValueDensity(NULL),fValueCT(NULL),fReadCalibration(false),
89 fMergedSlices(NULL),fDriverFile(
"Data.dat"),fCt2DensityFile(
"CT2Density.dat")
106 G4int returnvalue = 0;
size_t rflag = 0;
112 rflag = std::fread( buffer, 1, 128, dicom );
115 rflag = std::fread( buffer, 1, 4, dicom );
117 if(std::strncmp(
"DICM", buffer, 4) != 0) {
118 std::fseek(dicom, 0, SEEK_SET);
124 short elementLength2;
127 unsigned long elementLength4;
138 rflag = std::fread(buffer, 2, 1, dicom);
141 rflag = std::fread(buffer, 2, 1, dicom);
145 G4int tagDictionary = readGroupId*0x10000 + readElementId;
148 if(tagDictionary == 0x7FE00010) {
151 rflag = std::fread(buffer,2,1,dicom);
153 rflag = std::fread(buffer,4,1,dicom);
159 rflag = std::fread(buffer,2,1,dicom);
164 if((elementLength2 == 0x424f ||
165 elementLength2 == 0x574f ||
166 elementLength2 == 0x464f ||
167 elementLength2 == 0x5455 ||
168 elementLength2 == 0x5153 ||
169 elementLength2 == 0x4e55) &&
172 rflag = std::fread(buffer, 2, 1, dicom);
175 rflag = std::fread(buffer, 4, 1, dicom);
178 if(elementLength2 == 0x5153)
180 if(elementLength4 == 0xFFFFFFFF)
189 "Function read_defined_nested() failed!");
194 rflag = std::fread(data, elementLength4,1,dicom);
204 rflag = std::fread(buffer, 2, 1, dicom);
206 elementLength4 = elementLength2;
208 rflag = std::fread(data, elementLength4, 1, dicom);
215 if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
222 rflag = std::fread(buffer, 4, 1, dicom);
227 if(elementLength4 == 0xFFFFFFFF)
232 rflag = std::fread(data, elementLength4, 1, dicom);
239 data[elementLength4] =
'\0';
250 std::map<G4float,G4String>::const_iterator ite;
287 if (rflag)
return returnvalue;
295 if(tagDictionary == 0x00280010 ) {
297 std::printf(
"[0x00280010] Rows -> %i\n",
fRows);
299 }
else if(tagDictionary == 0x00280011 ) {
301 std::printf(
"[0x00280011] Columns -> %i\n",
fColumns);
303 }
else if(tagDictionary == 0x00280102 ) {
306 std::printf(
"[0x00280102] High bits -> %i\n",highBits);
308 }
else if(tagDictionary == 0x00280100 ) {
310 std::printf(
"[0x00280100] Bits allocated -> %i\n",
fBitAllocated);
312 }
else if(tagDictionary == 0x00280101 ) {
315 std::printf(
"[0x00280101] Bits stored -> %i\n",bitStored);
317 }
else if(tagDictionary == 0x00280106 ) {
319 std::printf(
"[0x00280106] Min. pixel value -> %i\n",
fMinPixelValue);
321 }
else if(tagDictionary == 0x00280107 ) {
323 std::printf(
"[0x00280107] Max. pixel value -> %i\n",
fMaxPixelValue);
325 }
else if(tagDictionary == 0x00281053) {
327 std::printf(
"[0x00281053] Rescale Slope -> %d\n",
fRescaleSlope);
329 }
else if(tagDictionary == 0x00281052 ) {
331 std::printf(
"[0x00281052] Rescale Intercept -> %d\n",
334 }
else if(tagDictionary == 0x00280103 ) {
337 std::printf(
"[0x00280103] Pixel Representation -> %i\n",
340 std::printf(
"### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
341 std::printf(
"DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
342 std::printf(
"ERROR !!!!!! -> \n");
345 }
else if(tagDictionary == 0x00080006 ) {
346 std::printf(
"[0x00080006] Modality -> %s\n", data);
348 }
else if(tagDictionary == 0x00080070 ) {
349 std::printf(
"[0x00080070] Manufacturer -> %s\n", data);
351 }
else if(tagDictionary == 0x00080080 ) {
352 std::printf(
"[0x00080080] Institution Name -> %s\n", data);
354 }
else if(tagDictionary == 0x00080081 ) {
355 std::printf(
"[0x00080081] Institution Address -> %s\n", data);
357 }
else if(tagDictionary == 0x00081040 ) {
358 std::printf(
"[0x00081040] Institution Department Name -> %s\n", data);
360 }
else if(tagDictionary == 0x00081090 ) {
361 std::printf(
"[0x00081090] Manufacturer's Model Name -> %s\n", data);
363 }
else if(tagDictionary == 0x00181000 ) {
364 std::printf(
"[0x00181000] Device Serial Number -> %s\n", data);
366 }
else if(tagDictionary == 0x00080008 ) {
367 std::printf(
"[0x00080008] Image Types -> %s\n", data);
369 }
else if(tagDictionary == 0x00283000 ) {
370 std::printf(
"[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
372 }
else if(tagDictionary == 0x00283002 ) {
373 std::printf(
"[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
375 }
else if(tagDictionary == 0x00283003 ) {
376 std::printf(
"[0x00283003] LUT Explanation LO 1 -> %s\n", data);
378 }
else if(tagDictionary == 0x00283004 ) {
379 std::printf(
"[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
381 }
else if(tagDictionary == 0x00283006 ) {
382 std::printf(
"[0x00283006] LUT Data US or SS -> %s\n", data);
384 }
else if(tagDictionary == 0x00283010 ) {
385 std::printf(
"[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
387 }
else if(tagDictionary == 0x00280120 ) {
388 std::printf(
"[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
390 }
else if(tagDictionary == 0x00280030 ) {
392 int iss = datas.find(
'\\');
394 fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() );
396 }
else if(tagDictionary == 0x00200037 ) {
397 std::printf(
"[0x00200037] Image Orientation (Phantom) -> %s\n", data);
399 }
else if(tagDictionary == 0x00200032 ) {
400 std::printf(
"[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
402 }
else if(tagDictionary == 0x00180050 ) {
404 std::printf(
"[0x00180050] Slice Thickness (mm) -> %f\n",
fSliceThickness);
406 }
else if(tagDictionary == 0x00201041 ) {
408 std::printf(
"[0x00201041] Slice Location -> %f\n",
fSliceLocation);
410 }
else if(tagDictionary == 0x00280004 ) {
412 std::printf(
"[0x00280004] Photometric Interpretation -> %s\n", data);
414 }
else if(tagDictionary == 0x00020010) {
415 if(strcmp(data,
"1.2.840.10008.1.2") == 0)
417 else if(strncmp(data,
"1.2.840.10008.1.2.2", 19) == 0)
421 std::printf(
"[0x00020010] Endian -> %s\n", data);
440 if(!dcmPZSH) {
return; }
466 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
467 mean +=
fTab[ww+sumy][xx+sumx];
514 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
515 mean +=
fTab[ww+sumy][xx+sumx];
537 foutG4DCM << density <<
" ";
538 if( xx%8 == 3 ) foutG4DCM <<
G4endl;
551 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
552 mean +=
fTab[ww+sumy][xx+sumx];
560 foutG4DCM << density <<
" ";
561 if( xx/fCompression%8 == 3 ) foutG4DCM <<
G4endl;
578 if( finData.eof() )
return;
581 for(
unsigned int ii = 0; ii < nMate; ii++ ){
582 finData >> mateName >> densityMax;
584 G4cout << ii <<
" ReadMaterialIndices " << mateName <<
" "
594 std::map<G4float,G4String>::reverse_iterator ite;
598 if( density >= (*ite).first ) {
615 G4int returnvalue = 0;
size_t rflag = 0;
627 std::printf(
"@@@ Error! Picture != 16 bits...\n");
628 std::printf(
"@@@ Error! Picture != 16 bits...\n");
629 std::printf(
"@@@ Error! Picture != 16 bits...\n");
631 unsigned char ch = 0;
636 rflag = std::fread( &ch, 1, 1, dicom);
648 rflag = std::fread(sbuff, 2, 1, dicom);
659 std::sprintf(nameProcessed,
"%s.g4dcmb",filename2);
660 fileOut = std::fopen(nameProcessed,
"w+b");
661 std::printf(
"### Writing of %s ###\n",nameProcessed);
664 rflag = std::fwrite(&nMate,
sizeof(
unsigned int), 1, fileOut);
666 std::map<G4float,G4String>::const_iterator ite;
669 for(
G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
673 const char* mateNameC = mateName.c_str();
674 rflag = std::fwrite(mateNameC,
sizeof(
char),40, fileOut);
679 unsigned int planesC = 1;
687 rflag = std::fwrite(&fRowsC,
sizeof(
unsigned int), 1, fileOut);
688 rflag = std::fwrite(&fColumnsC,
sizeof(
unsigned int), 1, fileOut);
689 rflag = std::fwrite(&planesC,
sizeof(
unsigned int), 1, fileOut);
691 rflag = std::fwrite(&pixelLocationXM,
sizeof(
G4float), 1, fileOut);
692 rflag = std::fwrite(&pixelLocationXP,
sizeof(
G4float), 1, fileOut);
693 rflag = std::fwrite(&pixelLocationYM,
sizeof(
G4float), 1, fileOut);
694 rflag = std::fwrite(&pixelLocationYP,
sizeof(
G4float), 1, fileOut);
695 rflag = std::fwrite(&fSliceLocationZM,
sizeof(
G4float), 1, fileOut);
696 rflag = std::fwrite(&fSliceLocationZP,
sizeof(
G4float), 1, fileOut);
699 std::printf(
"%8i %8i\n",fRows,
fColumns);
717 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
724 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
728 for(
int sumx = 0; sumx < compSize; sumx++) {
729 for(
int sumy = 0; sumy < compSize; sumy++) {
730 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
731 mean +=
fTab[ww+sumy][xx+sumx];
735 mean /= compSize*compSize;
740 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
753 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
760 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
764 for(
int sumx = 0; sumx < compSize; sumx++) {
765 for(
int sumy = 0; sumy < compSize; sumy++) {
766 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
767 mean +=
fTab[ww+sumy][xx+sumx];
771 mean /= compSize*compSize;
775 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
782 rflag = std::fclose(fileOut);
784 delete [] nameProcessed;
792 if (rflag)
return returnvalue;
816 "@@@ No value to transform pixels in density!");
847 density = fValueDensity[j] - ((fValueCT[j] - pixel)*deltaDensity/deltaCT );
853 std::printf(
"@@@ Error density = %f && Pixel = %i \
854 (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
855 deltaDensity/deltaCT);
866 char * oneLine =
new char[128];
868 if(!(checkData.is_open())) {
871 "\nDicomG4 needs Data.dat (or another driver file specified";
872 message +=
" in command line):\n";
873 message +=
"\tFirst line: number of image pixel for a voxel (G4Box)\n";
874 message +=
"\tSecond line: number of images (CT slices) to read\n";
875 message +=
"\tEach following line contains the name of a Dicom image";
876 message +=
" except for the .dcm extension";
886 checkData.getline(oneLine,100);
887 std::ifstream testExistence;
888 G4bool existAlready =
true;
890 checkData.getline(oneLine,100);
893 G4cout << fNFiles <<
" test file " << oneName <<
G4endl;
894 testExistence.open(oneName.
data());
895 if(!(testExistence.is_open())) {
896 existAlready =
false;
897 testExistence.clear();
898 testExistence.close();
900 testExistence.clear();
901 testExistence.close();
909 if( existAlready ==
false ) {
911 G4cout <<
"\nAll the necessary images were not found in processed form "
912 <<
", starting with .dcm images\n";
924 rflag = std::fscanf(lecturePref,
"%s",fCompressionc);
925 fCompression = atoi(fCompressionc);
926 rflag = std::fscanf(lecturePref,
"%s",maxc);
927 fNFiles = atoi(maxc);
932 rflag = std::fscanf(lecturePref,
"%s",inputFile);
933 std::sprintf(name,
"%s.dcm",inputFile);
936 std::printf(
"### Opening %s and reading :\n",name);
937 dicom = std::fopen(name,
"rb");
946 rflag = std::fclose(dicom);
948 rflag = std::fclose(lecturePref);
953 delete [] fCompressionc;
972 unsigned short item_GroupNumber;
973 unsigned short item_ElementNumber;
975 G4int items_array_length=0;
979 while(items_array_length < SQ_Length)
981 rflag = std::fread(buffer, 2, 1, nested);
984 rflag = std::fread(buffer, 2, 1, nested);
985 GetValue(buffer, item_ElementNumber);
987 rflag = std::fread(buffer, 4, 1, nested);
990 rflag = std::fread(buffer, item_Length, 1, nested);
992 items_array_length= items_array_length+8+item_Length;
997 if( SQ_Length>items_array_length )
1001 if (rflag)
return 1;
1009 unsigned short item_GroupNumber;
1010 unsigned short item_ElementNumber;
1011 unsigned int item_Length;
1017 rflag = std::fread(buffer, 2, 1, nested);
1018 GetValue(buffer, item_GroupNumber);
1020 rflag = std::fread(buffer, 2, 1, nested);
1021 GetValue(buffer, item_ElementNumber);
1023 rflag = std::fread(buffer, 4, 1, nested);
1026 if(item_Length!=0xffffffff)
1027 rflag = std::fread(buffer, item_Length, 1, nested);
1032 }
while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD
1044 unsigned short item_GroupNumber;
1045 unsigned short item_ElementNumber;
1046 G4int item_Length;
size_t rflag = 0;
1051 rflag = std::fread(buffer, 2, 1, nested);
1052 GetValue(buffer, item_GroupNumber);
1054 rflag = std::fread(buffer, 2, 1, nested);
1055 GetValue(buffer, item_ElementNumber);
1057 rflag = std::fread(buffer, 4, 1, nested);
1062 rflag = std::fread(buffer,item_Length,1,nested);
1065 while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D
1074 template <
class Type>
1077 #if BYTE_ORDER == BIG_ENDIAN
1079 #else // BYTE_ORDER == LITTLE_ENDIAN
1082 const int SIZE =
sizeof(_rval);
1084 for(
int i = 0; i < SIZE/2; i++) {
1086 _val[i] = _val[SIZE - 1 - i];
1087 _val[SIZE - 1 - i] = ctemp;
1090 _rval = *(Type *)_val;
static DicomHandler * fInstance
DicomHandler.cc :
void GetValue(char *, Type &)
void read_undefined_nested(FILE *)
void GetInformation(G4int &, char *)
const G4double w[NPOINTSGL]
Definition of the DicomPhantomZSliceMerged class.
static DicomHandler * Instance()
G4GLOB_DLL std::ostream G4cout
unsigned int GetMaterialIndex(G4float density)
G4int ReadFile(FILE *, char *)
void StoreData(std::ofstream &foutG4DCM)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
const char * data() const
G4int ReadData(FILE *, char *)
short fPixelRepresentation
Definition of the DicomHandler class.
G4float Pixel2density(G4int pixel)
std::map< G4float, G4String > fMaterialIndices
void AddZSlice(DicomPhantomZSliceHeader *val)
void ReadMaterialIndices(std::ifstream &finData)
G4int read_defined_nested(FILE *, G4int)
void read_undefined_item(FILE *)
DicomPhantomZSliceMerged * fMergedSlices