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), nbrequali(0),
88 valuedensity(NULL),valueCT(NULL),readCalibration(false),
89 mergedSlices(NULL),driverFile(
"Data.dat"),ct2densityFile(
"CT2Density.dat")
107 G4int returnvalue = 0;
size_t rflag = 0;
113 rflag = std::fread( buffer, 1, 128, dicom );
116 rflag = std::fread( buffer, 1, 4, dicom );
118 if(std::strncmp(
"DICM", buffer, 4) != 0) {
119 std::fseek(dicom, 0, SEEK_SET);
125 short elementLength2;
128 unsigned long elementLength4;
140 rflag = std::fread(buffer, 2, 1, dicom);
143 rflag = std::fread(buffer, 2, 1, dicom);
147 G4int tagDictionary = readGroupId*0x10000 + readElementId;
150 if(tagDictionary == 0x7FE00010) {
153 rflag = std::fread(buffer,2,1,dicom);
155 rflag = std::fread(buffer,4,1,dicom);
161 rflag = std::fread(buffer,2,1,dicom);
166 if((elementLength2 == 0x424f ||
167 elementLength2 == 0x574f ||
168 elementLength2 == 0x464f ||
169 elementLength2 == 0x5455 ||
170 elementLength2 == 0x5153 ||
171 elementLength2 == 0x4e55) &&
174 rflag = std::fread(buffer, 2, 1, dicom);
177 rflag = std::fread(buffer, 4, 1, dicom);
180 if(elementLength2 == 0x5153)
182 if(elementLength4 == 0xFFFFFFFF)
191 "Function read_defined_nested() failed!");
196 rflag = std::fread(data, elementLength4,1,dicom);
207 rflag = std::fread(buffer, 2, 1, dicom);
209 elementLength4 = elementLength2;
211 rflag = std::fread(data, elementLength4, 1, dicom);
218 if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
225 rflag = std::fread(buffer, 4, 1, dicom);
230 if(elementLength4 == 0xFFFFFFFF)
235 rflag = std::fread(data, elementLength4, 1, dicom);
242 data[elementLength4] =
'\0';
253 std::map<G4float,G4String>::const_iterator ite;
290 if (rflag)
return returnvalue;
298 if(tagDictionary == 0x00280010 ) {
300 std::printf(
"[0x00280010] Rows -> %i\n",
fRows);
302 }
else if(tagDictionary == 0x00280011 ) {
304 std::printf(
"[0x00280011] Columns -> %i\n",
fColumns);
306 }
else if(tagDictionary == 0x00280102 ) {
309 std::printf(
"[0x00280102] High bits -> %i\n",highBits);
311 }
else if(tagDictionary == 0x00280100 ) {
313 std::printf(
"[0x00280100] Bits allocated -> %i\n",
fBitAllocated);
315 }
else if(tagDictionary == 0x00280101 ) {
318 std::printf(
"[0x00280101] Bits stored -> %i\n",bitStored);
320 }
else if(tagDictionary == 0x00280106 ) {
322 std::printf(
"[0x00280106] Min. pixel value -> %i\n",
fMinPixelValue);
324 }
else if(tagDictionary == 0x00280107 ) {
326 std::printf(
"[0x00280107] Max. pixel value -> %i\n",
fMaxPixelValue);
328 }
else if(tagDictionary == 0x00281053) {
330 std::printf(
"[0x00281053] Rescale Slope -> %d\n",
fRescaleSlope);
332 }
else if(tagDictionary == 0x00281052 ) {
334 std::printf(
"[0x00281052] Rescale Intercept -> %d\n",
337 }
else if(tagDictionary == 0x00280103 ) {
340 std::printf(
"[0x00280103] Pixel Representation -> %i\n",
343 std::printf(
"### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
344 std::printf(
"DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
345 std::printf(
"ERROR !!!!!! -> \n");
348 }
else if(tagDictionary == 0x00080006 ) {
349 std::printf(
"[0x00080006] Modality -> %s\n", data);
351 }
else if(tagDictionary == 0x00080070 ) {
352 std::printf(
"[0x00080070] Manufacturer -> %s\n", data);
354 }
else if(tagDictionary == 0x00080080 ) {
355 std::printf(
"[0x00080080] Institution Name -> %s\n", data);
357 }
else if(tagDictionary == 0x00080081 ) {
358 std::printf(
"[0x00080081] Institution Address -> %s\n", data);
360 }
else if(tagDictionary == 0x00081040 ) {
361 std::printf(
"[0x00081040] Institution Department Name -> %s\n", data);
363 }
else if(tagDictionary == 0x00081090 ) {
364 std::printf(
"[0x00081090] Manufacturer's Model Name -> %s\n", data);
366 }
else if(tagDictionary == 0x00181000 ) {
367 std::printf(
"[0x00181000] Device Serial Number -> %s\n", data);
369 }
else if(tagDictionary == 0x00080008 ) {
370 std::printf(
"[0x00080008] Image Types -> %s\n", data);
372 }
else if(tagDictionary == 0x00283000 ) {
373 std::printf(
"[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
375 }
else if(tagDictionary == 0x00283002 ) {
376 std::printf(
"[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
378 }
else if(tagDictionary == 0x00283003 ) {
379 std::printf(
"[0x00283003] LUT Explanation LO 1 -> %s\n", data);
381 }
else if(tagDictionary == 0x00283004 ) {
382 std::printf(
"[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
384 }
else if(tagDictionary == 0x00283006 ) {
385 std::printf(
"[0x00283006] LUT Data US or SS -> %s\n", data);
387 }
else if(tagDictionary == 0x00283010 ) {
388 std::printf(
"[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
390 }
else if(tagDictionary == 0x00280120 ) {
391 std::printf(
"[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
393 }
else if(tagDictionary == 0x00280030 ) {
395 int iss = datas.find(
'\\');
397 fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() );
399 }
else if(tagDictionary == 0x00200037 ) {
400 std::printf(
"[0x00200037] Image Orientation (Phantom) -> %s\n", data);
402 }
else if(tagDictionary == 0x00200032 ) {
403 std::printf(
"[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
405 }
else if(tagDictionary == 0x00180050 ) {
407 std::printf(
"[0x00180050] Slice Thickness (mm) -> %f\n",
fSliceThickness);
409 }
else if(tagDictionary == 0x00201041 ) {
411 std::printf(
"[0x00201041] Slice Location -> %f\n",
fSliceLocation);
413 }
else if(tagDictionary == 0x00280004 ) {
415 std::printf(
"[0x00280004] Photometric Interpretation -> %s\n", data);
417 }
else if(tagDictionary == 0x00020010) {
418 if(strcmp(data,
"1.2.840.10008.1.2") == 0)
420 else if(strncmp(data,
"1.2.840.10008.1.2.2", 19) == 0)
424 std::printf(
"[0x00020010] Endian -> %s\n", data);
443 if(!dcmPZSH) {
return; }
469 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
470 mean +=
fTab[ww+sumy][xx+sumx];
517 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
518 mean +=
fTab[ww+sumy][xx+sumx];
540 foutG4DCM << density <<
" ";
541 if( xx%8 == 3 ) foutG4DCM <<
G4endl;
554 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
555 mean +=
fTab[ww+sumy][xx+sumx];
563 foutG4DCM << density <<
" ";
564 if( xx/fCompression%8 == 3 ) foutG4DCM <<
G4endl;
581 if( finData.eof() )
return;
584 for(
unsigned int ii = 0; ii < nMate; ii++ ){
585 finData >> mateName >> densityMax;
587 G4cout << ii <<
" ReadMaterialIndices " << mateName <<
" "
597 std::map<G4float,G4String>::reverse_iterator ite;
601 if( density >= (*ite).first ) {
618 G4int returnvalue = 0;
size_t rflag = 0;
630 std::printf(
"@@@ Error! Picture != 16 bits...\n");
631 std::printf(
"@@@ Error! Picture != 16 bits...\n");
632 std::printf(
"@@@ Error! Picture != 16 bits...\n");
634 unsigned char ch = 0;
639 rflag = std::fread( &ch, 1, 1, dicom);
651 rflag = std::fread(sbuff, 2, 1, dicom);
662 std::sprintf(nameProcessed,
"%s.g4dcmb",filename2);
663 fileOut = std::fopen(nameProcessed,
"w+b");
664 std::printf(
"### Writing of %s ###\n",nameProcessed);
667 rflag = std::fwrite(&nMate,
sizeof(
unsigned int), 1, fileOut);
669 std::map<G4float,G4String>::const_iterator ite;
672 for(
G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
676 const char* mateNameC = mateName.c_str();
677 rflag = std::fwrite(mateNameC,
sizeof(
char),40, fileOut);
682 unsigned int planesC = 1;
690 rflag = std::fwrite(&fRowsC,
sizeof(
unsigned int), 1, fileOut);
691 rflag = std::fwrite(&fColumnsC,
sizeof(
unsigned int), 1, fileOut);
692 rflag = std::fwrite(&planesC,
sizeof(
unsigned int), 1, fileOut);
694 rflag = std::fwrite(&pixelLocationXM,
sizeof(
G4float), 1, fileOut);
695 rflag = std::fwrite(&pixelLocationXP,
sizeof(
G4float), 1, fileOut);
696 rflag = std::fwrite(&pixelLocationYM,
sizeof(
G4float), 1, fileOut);
697 rflag = std::fwrite(&pixelLocationYP,
sizeof(
G4float), 1, fileOut);
698 rflag = std::fwrite(&fSliceLocationZM,
sizeof(
G4float), 1, fileOut);
699 rflag = std::fwrite(&fSliceLocationZP,
sizeof(
G4float), 1, fileOut);
702 std::printf(
"%8i %8i\n",fRows,
fColumns);
720 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
727 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
731 for(
int sumx = 0; sumx < compSize; sumx++) {
732 for(
int sumy = 0; sumy < compSize; sumy++) {
733 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
734 mean +=
fTab[ww+sumy][xx+sumx];
738 mean /= compSize*compSize;
743 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
756 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
763 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
767 for(
int sumx = 0; sumx < compSize; sumx++) {
768 for(
int sumy = 0; sumy < compSize; sumy++) {
769 if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow =
true;
770 mean +=
fTab[ww+sumy][xx+sumx];
774 mean /= compSize*compSize;
778 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
785 rflag = std::fclose(fileOut);
787 delete [] nameProcessed;
795 if (rflag)
return returnvalue;
820 "@@@ No value to transform pixels in density!");
851 density = valuedensity[j] - ((valueCT[j] - pixel)*deltaDensity/deltaCT );
857 std::printf(
"@@@ Error density = %f && Pixel = %i \
858 (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
859 deltaDensity/deltaCT);
870 char * oneLine =
new char[128];
872 if(!(checkData.is_open())) {
875 "\nDicomG4 needs Data.dat (or another driver file specified";
876 message +=
" in command line):\n";
877 message +=
"\tFirst line: number of image pixel for a voxel (G4Box)\n";
878 message +=
"\tSecond line: number of images (CT slices) to read\n";
879 message +=
"\tEach following line contains the name of a Dicom image";
880 message +=
" except for the .dcm extension";
890 checkData.getline(oneLine,100);
891 std::ifstream testExistence;
892 G4bool existAlready =
true;
894 checkData.getline(oneLine,100);
897 G4cout << fNFiles <<
" test file " << oneName <<
G4endl;
898 testExistence.open(oneName.
data());
899 if(!(testExistence.is_open())) {
900 existAlready =
false;
901 testExistence.clear();
902 testExistence.close();
904 testExistence.clear();
905 testExistence.close();
913 if( existAlready ==
false ) {
915 G4cout <<
"\nAll the necessary images were not found in processed form "
916 <<
", starting with .dcm images\n";
927 lecturePref = std::fopen(
driverFile.c_str(),
"r");
928 rflag = std::fscanf(lecturePref,
"%s",fCompressionc);
929 fCompression = atoi(fCompressionc);
930 rflag = std::fscanf(lecturePref,
"%s",maxc);
931 fNFiles = atoi(maxc);
936 rflag = std::fscanf(lecturePref,
"%s",inputFile);
937 std::sprintf(name,
"%s.dcm",inputFile);
940 std::printf(
"### Opening %s and reading :\n",name);
941 dicom = std::fopen(name,
"rb");
950 rflag = std::fclose(dicom);
952 rflag = std::fclose(lecturePref);
957 delete [] fCompressionc;
976 unsigned short item_GroupNumber;
977 unsigned short item_ElementNumber;
979 G4int items_array_length=0;
983 while(items_array_length < SQ_Length)
985 rflag = std::fread(buffer, 2, 1, nested);
988 rflag = std::fread(buffer, 2, 1, nested);
989 GetValue(buffer, item_ElementNumber);
991 rflag = std::fread(buffer, 4, 1, nested);
994 rflag = std::fread(buffer, item_Length, 1, nested);
996 items_array_length= items_array_length+8+item_Length;
1001 if( SQ_Length>items_array_length )
1005 if (rflag)
return 1;
1013 unsigned short item_GroupNumber;
1014 unsigned short item_ElementNumber;
1015 unsigned int item_Length;
1021 rflag = std::fread(buffer, 2, 1, nested);
1022 GetValue(buffer, item_GroupNumber);
1024 rflag = std::fread(buffer, 2, 1, nested);
1025 GetValue(buffer, item_ElementNumber);
1027 rflag = std::fread(buffer, 4, 1, nested);
1030 if(item_Length!=0xffffffff)
1031 rflag = std::fread(buffer, item_Length, 1, nested);
1036 }
while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD
1048 unsigned short item_GroupNumber;
1049 unsigned short item_ElementNumber;
1050 G4int item_Length;
size_t rflag = 0;
1055 rflag = std::fread(buffer, 2, 1, nested);
1056 GetValue(buffer, item_GroupNumber);
1058 rflag = std::fread(buffer, 2, 1, nested);
1059 GetValue(buffer, item_ElementNumber);
1061 rflag = std::fread(buffer, 4, 1, nested);
1066 rflag = std::fread(buffer,item_Length,1,nested);
1069 while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D
1078 template <
class Type>
1081 #if BYTE_ORDER == BIG_ENDIAN
1083 #else // BYTE_ORDER == LITTLE_ENDIAN
1086 const int SIZE =
sizeof(_rval);
1088 for(
int i = 0; i < SIZE/2; i++) {
1090 _val[i] = _val[SIZE - 1 - i];
1091 _val[SIZE - 1 - i] = ctemp;
1094 _rval = *(Type *)_val;
void GetValue(char *, Type &)
void read_undefined_nested(FILE *)
void GetInformation(G4int &, char *)
Definition of the DicomPhantomZSliceMerged class.
DicomPhantomZSliceMerged * mergedSlices
static DicomHandler * fgInstance
DicomHandler.cc :
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 *)