140 fNbrequali(0),fValueDensity(NULL),fValueCT(NULL),fReadCalibration(false),
149:DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
150 fCompression(0), fNFiles(0), fRows(0), fColumns(0),
151 fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
152 fPixelSpacingX(0.), fPixelSpacingY(0.),fSliceThickness(0.),
153 fSliceLocation(0.), fRescaleIntercept(0), fRescaleSlope(0),
154 fLittleEndian(true),fImplicitEndian(false),fPixelRepresentation(0),
155 fNbrequali(0),fValueDensity(NULL),fValueCT(NULL),
156 fReadCalibration(false),fMergedSlices(NULL),
157 fDriverFile(
"Data.dat"),fCt2DensityFile(
"CT2Density.dat")
171 G4cout <<
" ReadFile " << filename2 << G4endl;
173 G4int returnvalue = 0;
size_t rflag = 0;
179 rflag = std::fread( buffer, 1, 128, dicom );
182 rflag = std::fread( buffer, 1, 4, dicom );
184 if(std::strncmp(
"DICM", buffer, 4) != 0) {
185 std::fseek(dicom, 0, SEEK_SET);
191 short elementLength2;
194 unsigned long elementLength4;
205 rflag = std::fread(buffer, 2, 1, dicom);
208 rflag = std::fread(buffer, 2, 1, dicom);
212 G4int tagDictionary = readGroupId*0x10000 + readElementId;
215 if(tagDictionary == 0x7FE00010) {
219 rflag = std::fread(buffer,2,1,dicom);
221 rflag = std::fread(buffer,4,1,dicom);
227 rflag = std::fread(buffer,2,1,dicom);
232 if((elementLength2 == 0x424f ||
233 elementLength2 == 0x574f ||
234 elementLength2 == 0x464f ||
235 elementLength2 == 0x5455 ||
236 elementLength2 == 0x5153 ||
237 elementLength2 == 0x4e55) &&
240 rflag = std::fread(buffer, 2, 1, dicom);
243 rflag = std::fread(buffer, 4, 1, dicom);
246 if(elementLength2 == 0x5153)
248 if(elementLength4 == 0xFFFFFFFF)
254 G4Exception(
"DicomHandler::ReadFile",
257 "Function read_defined_nested() failed!");
262 rflag = std::fread(data, elementLength4,1,dicom);
272 rflag = std::fread(buffer, 2, 1, dicom);
274 elementLength4 = elementLength2;
276 rflag = std::fread(data, elementLength4, 1, dicom);
283 if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
284 G4Exception(
"DicomHandler::ReadFile",
290 rflag = std::fread(buffer, 4, 1, dicom);
295 if(elementLength4 == 0xFFFFFFFF)
300 rflag = std::fread(data, elementLength4, 1, dicom);
307 data[elementLength4] =
'\0';
355 if (rflag)
return returnvalue;
363 if(tagDictionary == 0x00280010 ) {
365 std::printf(
"[0x00280010] Rows -> %i\n",
fRows);
367 }
else if(tagDictionary == 0x00280011 ) {
369 std::printf(
"[0x00280011] Columns -> %i\n",
fColumns);
371 }
else if(tagDictionary == 0x00280102 ) {
374 std::printf(
"[0x00280102] High bits -> %i\n",highBits);
376 }
else if(tagDictionary == 0x00280100 ) {
378 std::printf(
"[0x00280100] Bits allocated -> %i\n",
fBitAllocated);
380 }
else if(tagDictionary == 0x00280101 ) {
383 std::printf(
"[0x00280101] Bits stored -> %i\n",bitStored);
385 }
else if(tagDictionary == 0x00280106 ) {
387 std::printf(
"[0x00280106] Min. pixel value -> %i\n",
fMinPixelValue);
389 }
else if(tagDictionary == 0x00280107 ) {
391 std::printf(
"[0x00280107] Max. pixel value -> %i\n",
fMaxPixelValue);
393 }
else if(tagDictionary == 0x00281053) {
395 std::printf(
"[0x00281053] Rescale Slope -> %d\n",
fRescaleSlope);
397 }
else if(tagDictionary == 0x00281052 ) {
399 std::printf(
"[0x00281052] Rescale Intercept -> %d\n",
402 }
else if(tagDictionary == 0x00280103 ) {
405 std::printf(
"[0x00280103] Pixel Representation -> %i\n",
408 std::printf(
"### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
409 std::printf(
"DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
410 std::printf(
"ERROR !!!!!! -> \n");
413 }
else if(tagDictionary == 0x00080006 ) {
414 std::printf(
"[0x00080006] Modality -> %s\n", data);
416 }
else if(tagDictionary == 0x00080070 ) {
417 std::printf(
"[0x00080070] Manufacturer -> %s\n", data);
419 }
else if(tagDictionary == 0x00080080 ) {
420 std::printf(
"[0x00080080] Institution Name -> %s\n", data);
422 }
else if(tagDictionary == 0x00080081 ) {
423 std::printf(
"[0x00080081] Institution Address -> %s\n", data);
425 }
else if(tagDictionary == 0x00081040 ) {
426 std::printf(
"[0x00081040] Institution Department Name -> %s\n", data);
428 }
else if(tagDictionary == 0x00081090 ) {
429 std::printf(
"[0x00081090] Manufacturer's Model Name -> %s\n", data);
431 }
else if(tagDictionary == 0x00181000 ) {
432 std::printf(
"[0x00181000] Device Serial Number -> %s\n", data);
434 }
else if(tagDictionary == 0x00080008 ) {
435 std::printf(
"[0x00080008] Image Types -> %s\n", data);
437 }
else if(tagDictionary == 0x00283000 ) {
438 std::printf(
"[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
440 }
else if(tagDictionary == 0x00283002 ) {
441 std::printf(
"[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
443 }
else if(tagDictionary == 0x00283003 ) {
444 std::printf(
"[0x00283003] LUT Explanation LO 1 -> %s\n", data);
446 }
else if(tagDictionary == 0x00283004 ) {
447 std::printf(
"[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
449 }
else if(tagDictionary == 0x00283006 ) {
450 std::printf(
"[0x00283006] LUT Data US or SS -> %s\n", data);
452 }
else if(tagDictionary == 0x00283010 ) {
453 std::printf(
"[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
455 }
else if(tagDictionary == 0x00280120 ) {
456 std::printf(
"[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
458 }
else if(tagDictionary == 0x00280030 ) {
460 G4int iss = G4int(datas.find(
'\\'));
462 fPixelSpacingY = atof( datas.substr(iss+1,datas.length()).c_str() );
464 }
else if(tagDictionary == 0x00200037 ) {
465 std::printf(
"[0x00200037] Image Orientation (Phantom) -> %s\n", data);
467 }
else if(tagDictionary == 0x00200032 ) {
468 std::printf(
"[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
470 }
else if(tagDictionary == 0x00180050 ) {
472 std::printf(
"[0x00180050] Slice Thickness (mm) -> %f\n",
fSliceThickness);
474 }
else if(tagDictionary == 0x00201041 ) {
476 std::printf(
"[0x00201041] Slice Location -> %f\n",
fSliceLocation);
478 }
else if(tagDictionary == 0x00280004 ) {
480 std::printf(
"[0x00280004] Photometric Interpretation -> %s\n", data);
482 }
else if(tagDictionary == 0x00020010) {
483 if(strcmp(data,
"1.2.840.10008.1.2") == 0)
485 else if(strncmp(data,
"1.2.840.10008.1.2.2", 19) == 0)
489 std::printf(
"[0x00020010] Endian -> %s\n", data);
681 G4int returnvalue = 0;
size_t rflag = 0;
686 for ( G4int i = 0; i <
fRows; ++i )
693 std::printf(
"@@@ Error! Picture != 16 bits...\n");
694 std::printf(
"@@@ Error! Picture != 16 bits...\n");
695 std::printf(
"@@@ Error! Picture != 16 bits...\n");
697 unsigned char ch = 0;
699 for(G4int j = 0; j <
fRows; ++j) {
700 for(G4int i = 0; i <
fColumns; ++i) {
701 rflag = std::fread( &ch, 1, 1, dicom);
710 for( G4int j = 0; j <
fRows; ++j) {
711 for( G4int i = 0; i <
fColumns; ++i) {
712 rflag = std::fread(sbuff, 2, 1, dicom);
721 FILE* fileOut = std::fopen(nameProcessed.c_str(),
"w+b");
723 G4cout <<
"### Writing of "<< nameProcessed <<
" ###\n";
726 rflag = std::fwrite(&nMate,
sizeof(
unsigned int), 1, fileOut);
732 for( std::size_t ii = (*ite).second.length(); ii < 40; ++ii )
737 const char* mateNameC = mateName.c_str();
738 rflag = std::fwrite(mateNameC,
sizeof(
char),40, fileOut);
743 unsigned int planesC = 1;
751 rflag = std::fwrite(&fRowsC,
sizeof(
unsigned int), 1, fileOut);
752 rflag = std::fwrite(&fColumnsC,
sizeof(
unsigned int), 1, fileOut);
753 rflag = std::fwrite(&planesC,
sizeof(
unsigned int), 1, fileOut);
755 rflag = std::fwrite(&pixelLocationXM,
sizeof(G4float), 1, fileOut);
756 rflag = std::fwrite(&pixelLocationXP,
sizeof(G4float), 1, fileOut);
757 rflag = std::fwrite(&pixelLocationYM,
sizeof(G4float), 1, fileOut);
758 rflag = std::fwrite(&pixelLocationYP,
sizeof(G4float), 1, fileOut);
759 rflag = std::fwrite(&fSliceLocationZM,
sizeof(G4float), 1, fileOut);
760 rflag = std::fwrite(&fSliceLocationZP,
sizeof(G4float), 1, fileOut);
772 G4bool overflow =
false;
776 for( G4int ww = 0; ww <
fRows; ++ww) {
777 for( G4int xx = 0; xx <
fColumns; ++xx) {
781 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
788 for(G4int ww = 0; ww <
fRows ;ww += compSize ) {
789 for(G4int xx = 0; xx <
fColumns ;xx +=compSize ) {
792 for(G4int sumx = 0; sumx < compSize; ++sumx) {
793 for(G4int sumy = 0; sumy < compSize; ++sumy) {
795 mean +=
fTab[ww+sumy][xx+sumx];
799 mean /= compSize*compSize;
804 rflag = std::fwrite(&mateID,
sizeof(
unsigned int), 1, fileOut);
813 for( G4int ww = 0; ww <
fRows; ++ww) {
814 for( G4int xx = 0; xx <
fColumns; ++xx) {
817 rflag = std::fwrite(&density,
sizeof(G4float), 1, fileOut);
824 for(G4int ww = 0; ww <
fRows ;ww += compSize ) {
825 for(G4int xx = 0; xx <
fColumns ;xx +=compSize ) {
828 for(G4int sumx = 0; sumx < compSize; ++sumx) {
829 for(G4int sumy = 0; sumy < compSize; ++sumy) {
831 mean +=
fTab[ww+sumy][xx+sumx];
835 mean /= compSize*compSize;
839 rflag = std::fwrite(&density,
sizeof(G4float), 1, fileOut);
846 rflag = std::fclose(fileOut);
848 delete [] nameProcessed;
856 if (rflag)
return returnvalue;
965 char * oneLine =
new char[128];
967 if(!(checkData.is_open())) {
970 "\nDicomG4 needs Data.dat (or another driver file specified";
971 message +=
" in command line):\n";
972 message +=
"\tFirst line: number of image pixel for a voxel (G4Box)\n";
973 message +=
"\tSecond line: number of images (CT slices) to read\n";
974 message +=
"\tEach following line contains the name of a Dicom image";
975 message +=
" except for the .dcm extension";
976 G4Exception(
"DicomHandler::ReadFile",
985 checkData.getline(oneLine,100);
986 std::ifstream testExistence;
987 G4bool existAlready =
true;
988 for(G4int rep = 0; rep <
fNFiles; ++rep) {
989 checkData.getline(oneLine,100);
992 G4cout <<
fNFiles <<
" test file " << oneName << G4endl;
993 testExistence.open(oneName.data());
994 if(!(testExistence.is_open())) {
995 existAlready =
false;
996 testExistence.clear();
997 testExistence.close();
999 testExistence.clear();
1000 testExistence.close();
1008 if( existAlready ==
false ) {
1010 G4cout <<
"\nAll the necessary images were not found in processed form "
1011 <<
", starting with .dcm images\n";
1020 lecturePref = std::fopen(
fDriverFile.c_str(),
"r");
1022 rflag = std::fscanf(lecturePref,
"%s",fCompressionc);
1024 rflag = std::fscanf(lecturePref,
"%s",maxc);
1026 G4cout <<
" fNFiles " <<
fNFiles << G4endl;
1030#ifdef DICOM_USE_HEAD
1031 for( G4int i = 1; i <=
fNFiles; ++i )
1033 rflag = std::fscanf(lecturePref,
"%s",inputFile);
1041 G4cout <<
"### Opening " << name2 <<
" and reading :\n";
1042 dicom = std::fopen(name2.c_str(),
"rb");
1049 G4cout <<
"\nError opening file : " << name2 << G4endl;
1051 rflag = std::fclose(dicom);
1055 for( G4int i = 1; i <=
fNFiles; ++i )
1057 rflag = std::fscanf(lecturePref,
"%s",inputFile);
1064 G4cout <<
"### Opening " << name <<
" and reading :\n";
1065 dicom = std::fopen(name.c_str(),
"rb");
1072 G4cout <<
"\nError opening file : " << name << G4endl;
1074 rflag = std::fclose(dicom);
1078 rflag = std::fclose(lecturePref);
1083 delete [] fCompressionc;
1085 delete [] inputFile;
1101 unsigned short item_GroupNumber;
1102 unsigned short item_ElementNumber;
1104 G4int items_array_length=0;
1108 while(items_array_length < SQ_Length)
1110 rflag = std::fread(buffer, 2, 1, nested);
1111 GetValue(buffer, item_GroupNumber);
1113 rflag = std::fread(buffer, 2, 1, nested);
1114 GetValue(buffer, item_ElementNumber);
1116 rflag = std::fread(buffer, 4, 1, nested);
1119 rflag = std::fread(buffer, item_Length, 1, nested);
1121 items_array_length= items_array_length+8+item_Length;
1126 if( SQ_Length>items_array_length )
1130 if (rflag)
return 1;
1138 unsigned short item_GroupNumber;
1139 unsigned short item_ElementNumber;
1140 unsigned int item_Length;
1146 rflag = std::fread(buffer, 2, 1, nested);
1147 GetValue(buffer, item_GroupNumber);
1149 rflag = std::fread(buffer, 2, 1, nested);
1150 GetValue(buffer, item_ElementNumber);
1152 rflag = std::fread(buffer, 4, 1, nested);
1155 if(item_Length!=0xffffffff)
1156 rflag = std::fread(buffer, item_Length, 1, nested);
1161 }
while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD
1173 unsigned short item_GroupNumber;
1174 unsigned short item_ElementNumber;
1175 G4int item_Length;
size_t rflag = 0;
1180 rflag = std::fread(buffer, 2, 1, nested);
1181 GetValue(buffer, item_GroupNumber);
1183 rflag = std::fread(buffer, 2, 1, nested);
1184 GetValue(buffer, item_ElementNumber);
1186 rflag = std::fread(buffer, 4, 1, nested);
1191 rflag = std::fread(buffer,item_Length,1,nested);
1194 while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D