79 :DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
80 fCompression(0),fNFiles(0), fRows(0), fColumns(0),
81 fBitAllocated(0), fMaxPixelValue(0),
82 fMinPixelValue(0),fPixelSpacingX(0.),
83 fPixelSpacingY(0.),fSliceThickness(0.),
84 fSliceLocation(0.),fRescaleIntercept(0),
85 fRescaleSlope(0),fLittleEndian(true),
86 fImplicitEndian(false),fPixelRepresentation(0),
87 fNbrequali(0),fValueDensity(NULL),fValueCT(NULL),fReadCalibration(false),
88 fMergedSlices(NULL),fCt2DensityFile(
"null.dat")
91 G4String path = getenv(
"DICOM_PATH");
97 :DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
98 fCompression(0), fNFiles(0), fRows(0), fColumns(0),
99 fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
100 fPixelSpacingX(0.), fPixelSpacingY(0.),fSliceThickness(0.),
101 fSliceLocation(0.), fRescaleIntercept(0), fRescaleSlope(0),
102 fLittleEndian(true),fImplicitEndian(false),fPixelRepresentation(0),
103 fNbrequali(0),fValueDensity(NULL),fValueCT(NULL),
104 fReadCalibration(false),fMergedSlices(NULL),
105 fDriverFile(
"Data.dat"),fCt2DensityFile(
"CT2Density.dat")
121 G4int returnvalue = 0;
size_t rflag = 0;
127 rflag = std::fread( buffer, 1, 128, dicom );
130 rflag = std::fread( buffer, 1, 4, dicom );
132 if(std::strncmp(
"DICM", buffer, 4) != 0) {
133 std::fseek(dicom, 0, SEEK_SET);
139 short elementLength2;
142 unsigned long elementLength4;
153 rflag = std::fread(buffer, 2, 1, dicom);
156 rflag = std::fread(buffer, 2, 1, dicom);
160 G4int tagDictionary = readGroupId*0x10000 + readElementId;
163 if(tagDictionary == 0x7FE00010) {
166 rflag = std::fread(buffer,2,1,dicom);
168 rflag = std::fread(buffer,4,1,dicom);
174 rflag = std::fread(buffer,2,1,dicom);
179 if((elementLength2 == 0x424f ||
180 elementLength2 == 0x574f ||
181 elementLength2 == 0x464f ||
182 elementLength2 == 0x5455 ||
183 elementLength2 == 0x5153 ||
184 elementLength2 == 0x4e55) &&
187 rflag = std::fread(buffer, 2, 1, dicom);
190 rflag = std::fread(buffer, 4, 1, dicom);
193 if(elementLength2 == 0x5153)
195 if(elementLength4 == 0xFFFFFFFF)
204 "Function read_defined_nested() failed!");
209 rflag = std::fread(data, elementLength4,1,dicom);
219 rflag = std::fread(buffer, 2, 1, dicom);
221 elementLength4 = elementLength2;
223 rflag = std::fread(data, elementLength4, 1, dicom);
230 if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
237 rflag = std::fread(buffer, 4, 1, dicom);
242 if(elementLength4 == 0xFFFFFFFF)
247 rflag = std::fread(data, elementLength4, 1, dicom);
254 data[elementLength4] =
'\0';
264 std::map<G4float,G4String>::const_iterator ite;
301 if (rflag)
return returnvalue;
309 if(tagDictionary == 0x00280010 ) {
313 }
else if(tagDictionary == 0x00280011 ) {
317 }
else if(tagDictionary == 0x00280102 ) {
320 std::printf(
"[0x00280102] High bits -> %i\n",highBits);
322 }
else if(tagDictionary == 0x00280100 ) {
326 }
else if(tagDictionary == 0x00280101 ) {
329 std::printf(
"[0x00280101] Bits stored -> %i\n",bitStored);
331 }
else if(tagDictionary == 0x00280106 ) {
335 }
else if(tagDictionary == 0x00280107 ) {
339 }
else if(tagDictionary == 0x00281053) {
343 }
else if(tagDictionary == 0x00281052 ) {
345 std::printf(
"[0x00281052] Rescale Intercept -> %d\n",
348 }
else if(tagDictionary == 0x00280103 ) {
351 std::printf(
"[0x00280103] Pixel Representation -> %i\n",
354 std::printf(
"### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
355 std::printf(
"DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
359 }
else if(tagDictionary == 0x00080006 ) {
360 std::printf(
"[0x00080006] Modality -> %s\n", data);
362 }
else if(tagDictionary == 0x00080070 ) {
363 std::printf(
"[0x00080070] Manufacturer -> %s\n", data);
365 }
else if(tagDictionary == 0x00080080 ) {
366 std::printf(
"[0x00080080] Institution Name -> %s\n", data);
368 }
else if(tagDictionary == 0x00080081 ) {
369 std::printf(
"[0x00080081] Institution Address -> %s\n", data);
371 }
else if(tagDictionary == 0x00081040 ) {
372 std::printf(
"[0x00081040] Institution Department Name -> %s\n", data);
374 }
else if(tagDictionary == 0x00081090 ) {
375 std::printf(
"[0x00081090] Manufacturer's Model Name -> %s\n", data);
377 }
else if(tagDictionary == 0x00181000 ) {
378 std::printf(
"[0x00181000] Device Serial Number -> %s\n", data);
380 }
else if(tagDictionary == 0x00080008 ) {
381 std::printf(
"[0x00080008] Image Types -> %s\n", data);
383 }
else if(tagDictionary == 0x00283000 ) {
384 std::printf(
"[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
386 }
else if(tagDictionary == 0x00283002 ) {
387 std::printf(
"[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
389 }
else if(tagDictionary == 0x00283003 ) {
390 std::printf(
"[0x00283003] LUT Explanation LO 1 -> %s\n", data);
392 }
else if(tagDictionary == 0x00283004 ) {
393 std::printf(
"[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
395 }
else if(tagDictionary == 0x00283006 ) {
396 std::printf(
"[0x00283006] LUT Data US or SS -> %s\n", data);
398 }
else if(tagDictionary == 0x00283010 ) {
399 std::printf(
"[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
401 }
else if(tagDictionary == 0x00280120 ) {
402 std::printf(
"[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
404 }
else if(tagDictionary == 0x00280030 ) {
406 int iss = datas.find(
'\\');
408 fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() );
410 }
else if(tagDictionary == 0x00200037 ) {
411 std::printf(
"[0x00200037] Image Orientation (Phantom) -> %s\n", data);
413 }
else if(tagDictionary == 0x00200032 ) {
414 std::printf(
"[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
416 }
else if(tagDictionary == 0x00180050 ) {
420 }
else if(tagDictionary == 0x00201041 ) {
424 }
else if(tagDictionary == 0x00280004 ) {
426 std::printf(
"[0x00280004] Photometric Interpretation -> %s\n", data);
428 }
else if(tagDictionary == 0x00020010) {
429 if(strcmp(data,
"1.2.840.10008.1.2") == 0)
431 else if(strncmp(data,
"1.2.840.10008.1.2.2", 19) == 0)
454 if(!dcmPZSH) {
return; }
480 if(ww+sumy >= fRows ||
xx+sumx >= fColumns) overflow =
true;
481 mean +=
fTab[ww+sumy][
xx+sumx];
528 if(ww+sumy >= fRows ||
xx+sumx >= fColumns) overflow =
true;
529 mean +=
fTab[ww+sumy][
xx+sumx];
551 foutG4DCM << density <<
" ";
552 if(
xx%8 == 3 ) foutG4DCM <<
G4endl;
565 if(ww+sumy >= fRows ||
xx+sumx >= fColumns) overflow =
true;
566 mean +=
fTab[ww+sumy][
xx+sumx];
574 foutG4DCM << density <<
" ";
575 if(
xx/fCompression%8 == 3 ) foutG4DCM <<
G4endl;
592 if( finData.eof() )
return;
595 for(
unsigned int ii = 0; ii < nMate; ii++ ){
596 finData >> mateName >> densityMax;
608 std::map<G4float,G4String>::reverse_iterator ite;
613 if( density >= (*ite).first ) {
630 G4int returnvalue = 0;
size_t rflag = 0;
646 unsigned char ch = 0;
651 rflag = std::fread( &ch, 1, 1, dicom);
663 rflag = std::fread(sbuff, 2, 1, dicom);
674 std::sprintf(nameProcessed,
"%s.g4dcmb",filename2);
675 fileOut = std::fopen(nameProcessed,
"w+b");
676 std::printf(
"### Writing of %s ###\n",nameProcessed);
679 rflag = std::fwrite(&nMate,
sizeof(
unsigned int), 1, fileOut);
681 std::map<G4float,G4String>::const_iterator ite;
684 for(
G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
688 const char* mateNameC = mateName.c_str();
689 rflag = std::fwrite(mateNameC,
sizeof(
char),40, fileOut);
694 unsigned int planesC = 1;
702 rflag = std::fwrite(&fRowsC,
sizeof(
unsigned int), 1, fileOut);
703 rflag = std::fwrite(&fColumnsC,
sizeof(
unsigned int), 1, fileOut);
704 rflag = std::fwrite(&planesC,
sizeof(
unsigned int), 1, fileOut);
706 rflag = std::fwrite(&pixelLocationXM,
sizeof(
G4float), 1, fileOut);
707 rflag = std::fwrite(&pixelLocationXP,
sizeof(
G4float), 1, fileOut);
708 rflag = std::fwrite(&pixelLocationYM,
sizeof(
G4float), 1, fileOut);
709 rflag = std::fwrite(&pixelLocationYP,
sizeof(
G4float), 1, fileOut);
710 rflag = std::fwrite(&fSliceLocationZM,
sizeof(
G4float), 1, fileOut);
711 rflag = std::fwrite(&fSliceLocationZP,
sizeof(
G4float), 1, fileOut);
732 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
739 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
743 for(
int sumx = 0; sumx < compSize; sumx++) {
744 for(
int sumy = 0; sumy < compSize; sumy++) {
745 if(ww+sumy >= fRows ||
xx+sumx >= fColumns) overflow =
true;
746 mean +=
fTab[ww+sumy][
xx+sumx];
750 mean /= compSize*compSize;
755 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
768 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
775 for(
G4int ww = 0; ww <
fRows ;ww += compSize ) {
779 for(
int sumx = 0; sumx < compSize; sumx++) {
780 for(
int sumy = 0; sumy < compSize; sumy++) {
781 if(ww+sumy >= fRows ||
xx+sumx >= fColumns) overflow =
true;
782 mean +=
fTab[ww+sumy][
xx+sumx];
786 mean /= compSize*compSize;
790 rflag = std::fwrite(&density,
sizeof(
G4float), 1, fileOut);
799 delete [] nameProcessed;
807 if (rflag)
return returnvalue;
815 #ifdef DICOM_USE_HEAD
820 G4cout <<
"No calibration curve for the DICOM HEAD needed!" <<
G4endl;
839 "@@@ No value to transform pixels in density!");
852 #ifdef DICOM_USE_HEAD
858 if (pixel == -998.) density = 0.001290;
860 else if ( pixel == 24.) density = 1.00;
862 else if ( pixel == 52.) density = 1.03;
864 else if ( pixel == 92.) density = 1.10;
866 else if ( pixel == 197.) density = 1.18;
868 else if ( pixel == 923.) density = 1.85;
870 else if ( pixel == 1280.) density = 2.14;
872 else if ( pixel == 2310.) density = 2.89;
896 density = fValueDensity[j]-((fValueCT[j] - pixel)*deltaDensity/deltaCT );
902 std::printf(
"@@@ Error density = %f && Pixel = %i \
903 (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
904 deltaDensity/deltaCT);
915 char * oneLine =
new char[128];
917 if(!(checkData.is_open())) {
920 "\nDicomG4 needs Data.dat (or another driver file specified";
921 message +=
" in command line):\n";
922 message +=
"\tFirst line: number of image pixel for a voxel (G4Box)\n";
923 message +=
"\tSecond line: number of images (CT slices) to read\n";
924 message +=
"\tEach following line contains the name of a Dicom image";
925 message +=
" except for the .dcm extension";
935 checkData.getline(oneLine,100);
936 std::ifstream testExistence;
937 G4bool existAlready =
true;
939 checkData.getline(oneLine,100);
942 G4cout << fNFiles <<
" test file " << oneName <<
G4endl;
943 testExistence.open(oneName.
data());
944 if(!(testExistence.is_open())) {
945 existAlready =
false;
946 testExistence.clear();
947 testExistence.close();
949 testExistence.clear();
950 testExistence.close();
958 if( existAlready ==
false ) {
960 G4cout <<
"\nAll the necessary images were not found in processed form "
961 <<
", starting with .dcm images\n";
973 rflag = std::fscanf(lecturePref,
"%s",fCompressionc);
974 fCompression = atoi(fCompressionc);
975 rflag = std::fscanf(lecturePref,
"%s",maxc);
976 fNFiles = atoi(maxc);
981 #ifdef DICOM_USE_HEAD
984 rflag = std::fscanf(lecturePref,
"%s",inputFile);
988 std::sprintf(name,
"%s.dcm",inputFile);
992 strcpy(name2, path.c_str());
995 std::printf(
"### Opening %s and reading :\n",name2);
996 dicom = std::fopen(name2,
"rb");
1003 G4cout <<
"\nError opening file : " << name2 <<
G4endl;
1011 rflag = std::fscanf(lecturePref,
"%s",inputFile);
1013 std::sprintf(name,
"%s.dcm",inputFile);
1018 std::printf(
"### Opening %s and reading :\n",name);
1019 dicom = std::fopen(name,
"rb");
1037 delete [] fCompressionc;
1040 delete [] inputFile;
1056 unsigned short item_GroupNumber;
1057 unsigned short item_ElementNumber;
1059 G4int items_array_length=0;
1063 while(items_array_length < SQ_Length)
1065 rflag = std::fread(buffer, 2, 1, nested);
1066 GetValue(buffer, item_GroupNumber);
1068 rflag = std::fread(buffer, 2, 1, nested);
1069 GetValue(buffer, item_ElementNumber);
1071 rflag = std::fread(buffer, 4, 1, nested);
1074 rflag = std::fread(buffer, item_Length, 1, nested);
1076 items_array_length= items_array_length+8+item_Length;
1081 if( SQ_Length>items_array_length )
1085 if (rflag)
return 1;
1093 unsigned short item_GroupNumber;
1094 unsigned short item_ElementNumber;
1095 unsigned int item_Length;
1101 rflag = std::fread(buffer, 2, 1, nested);
1102 GetValue(buffer, item_GroupNumber);
1104 rflag = std::fread(buffer, 2, 1, nested);
1105 GetValue(buffer, item_ElementNumber);
1107 rflag = std::fread(buffer, 4, 1, nested);
1110 if(item_Length!=0xffffffff)
1111 rflag = std::fread(buffer, item_Length, 1, nested);
1116 }
while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD
1128 unsigned short item_GroupNumber;
1129 unsigned short item_ElementNumber;
1130 G4int item_Length;
size_t rflag = 0;
1135 rflag = std::fread(buffer, 2, 1, nested);
1136 GetValue(buffer, item_GroupNumber);
1138 rflag = std::fread(buffer, 2, 1, nested);
1139 GetValue(buffer, item_ElementNumber);
1141 rflag = std::fread(buffer, 4, 1, nested);
1146 rflag = std::fread(buffer,item_Length,1,nested);
1149 while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D
1158 template <
class Type>
1161 #if BYTE_ORDER == BIG_ENDIAN
1163 #else // BYTE_ORDER == LITTLE_ENDIAN
1166 const int SIZE =
sizeof(_rval);
1168 for(
int i = 0; i < SIZE/2; i++) {
1170 _val[i] = _val[SIZE - 1 - i];
1171 _val[SIZE - 1 - i] = ctemp;
1174 _rval = *(Type *)_val;
Definition of the DicomHandler class.
void AddZSlice(DicomPhantomZSliceHeader *val)
G4int ReadData(FILE *, char *)
Definition of the DicomPhantomZSliceMerged class.
void message(RunManager *runmanager)
std::map< G4float, G4String > fMaterialIndices
void read_undefined_item(FILE *)
void read_undefined_nested(FILE *)
const XML_Char const XML_Char * data
short fPixelRepresentation
const char * data() const
printf("%d Experimental points found\n", nlines)
G4float Pixel2density(G4int pixel)
G4int ReadFile(FILE *, char *)
void GetValue(char *, Type &)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static DicomHandler * Instance()
static DicomHandler * fInstance
void StoreData(std::ofstream &foutG4DCM)
G4GLOB_DLL std::ostream G4cout
unsigned int GetMaterialIndex(G4float density)
G4int read_defined_nested(FILE *, G4int)
void GetInformation(G4int &, char *)
DicomPhantomZSliceMerged * fMergedSlices
void ReadMaterialIndices(std::ifstream &finData)