File indexing completed on 2025-04-10 08:06:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 #include "DicomFileMgr.hh"
0027
0028 #include "DicomFileCT.hh"
0029 #include "DicomFilePET.hh"
0030 #include "DicomFilePlan.hh"
0031 #include "DicomFileStructure.hh"
0032 #include "dcmtk/dcmdata/dcdeftag.h"
0033
0034 #include "G4UIcommand.hh"
0035 #include "G4tgrFileIn.hh"
0036
0037 DicomFileMgr* DicomFileMgr::theInstance = 0;
0038 int DicomFileMgr::verbose = 1;
0039
0040
0041 DicomFileMgr* DicomFileMgr::GetInstance()
0042 {
0043 if (!theInstance) {
0044 theInstance = new DicomFileMgr;
0045 }
0046 return theInstance;
0047 }
0048
0049
0050 DicomFileMgr::DicomFileMgr()
0051 {
0052 fCompression = 1.;
0053 theCTFileAll = 0;
0054 theStructureNCheck = 4;
0055 theStructureNMaxROI = 100;
0056 }
0057
0058
0059 void DicomFileMgr::Convert(G4String fileName)
0060 {
0061 G4tgrFileIn fin = G4tgrFileIn::GetInstance(fileName);
0062 std::vector<G4String> wl;
0063
0064 theFileOutName = "test.g4dcm";
0065 int ii;
0066 for (ii = 0;; ii++) {
0067 if (!fin.GetWordsInLine(wl)) break;
0068 if (wl[0] == ":COMPRESSION") {
0069 CheckNColumns(wl, 2);
0070 SetCompression(wl[1]);
0071 }
0072 else if (wl[0] == ":FILE") {
0073 CheckNColumns(wl, 2);
0074 G4cout << "@@@@@@@ Reading FILE: " << wl[1] << G4endl;
0075 AddFile(wl[1]);
0076 }
0077 else if (wl[0] == ":FILE_OUT") {
0078 CheckNColumns(wl, 2);
0079 theFileOutName = wl[1];
0080 }
0081 else if (wl[0] == ":MATE_DENS") {
0082 CheckNColumns(wl, 3);
0083 AddMaterialDensity(wl);
0084 }
0085 else if (wl[0] == ":MATE") {
0086 CheckNColumns(wl, 3);
0087 AddMaterial(wl);
0088 }
0089 else if (wl[0] == ":CT2D") {
0090 CheckNColumns(wl, 3);
0091 AddCT2Density(wl);
0092 }
0093 else {
0094 G4Exception("DICOM2G4", "Wrong argument", FatalErrorInArgument,
0095 G4String("UNKNOWN TAG IN FILE " + wl[0]).c_str());
0096 }
0097 }
0098
0099
0100 ProcessFiles();
0101 }
0102
0103
0104 void DicomFileMgr::CheckNColumns(std::vector<G4String> wl, size_t vsizeTh)
0105 {
0106 if (wl.size() != vsizeTh) {
0107 G4cerr << " Reading line " << G4endl;
0108 for (size_t ii = 0; ii < wl.size(); ii++) {
0109 G4cerr << wl[ii] << " ";
0110 }
0111 G4cerr << G4endl;
0112 G4Exception("DICOM2G4", "D2G0010", FatalErrorInArgument,
0113 ("Wrong number of columns in line " + std::to_string(wl.size()) + " <> "
0114 + std::to_string(vsizeTh))
0115 .c_str());
0116 }
0117 }
0118
0119
0120 void DicomFileMgr::SetCompression(G4String fComp)
0121 {
0122 fCompression = G4UIcommand::ConvertToDouble(fComp);
0123 }
0124
0125
0126 void DicomFileMgr::AddFile(G4String fileName)
0127 {
0128 DcmFileFormat dfile;
0129 if (!(dfile.loadFile(fileName.c_str())).good()) {
0130 G4Exception("DicomHandler::ReadFile", "dfile.loadFile", FatalErrorInArgument,
0131 ("Error reading file " + fileName).c_str());
0132 }
0133 DcmDataset* dset = dfile.getDataset();
0134
0135 OFString dModality;
0136 if (!dset->findAndGetOFString(DCM_Modality, dModality).good()) {
0137 G4Exception("DicomHandler::ReadData ", "", FatalException, " Have not read Modality");
0138 }
0139
0140 if (dModality == "CT") {
0141 DicomFileCT* df = new DicomFileCT(dset);
0142 df->ReadData();
0143 df->SetFileName(fileName);
0144
0145 theCTFiles[df->GetMaxZ()] = df;
0146 }
0147 else if (dModality == "RTSTRUCT") {
0148 DicomFileStructure* df = new DicomFileStructure(dset);
0149 df->ReadData();
0150 df->SetFileName(fileName);
0151
0152 theStructFiles.push_back(df);
0153 }
0154 else if (dModality == "RTPLAN") {
0155 DicomFilePlan* df = new DicomFilePlan(dset);
0156 df->ReadData();
0157 df->SetFileName(fileName);
0158
0159 thePlanFiles.push_back(df);
0160 }
0161 else if (dModality == "PT") {
0162 DicomFilePET* df = new DicomFilePET(dset);
0163 df->ReadData();
0164 df->SetFileName(fileName);
0165
0166 thePETFiles[df->GetMaxZ()] = df;
0167
0168 }
0169 else {
0170 G4Exception(
0171 "DicomFileMgr::AddFIle()", "DFM001", FatalErrorInArgument,
0172 (G4String("File is not of type CT or RTSTRUCT or RTPLAN, but: ") + dModality).c_str());
0173 }
0174 }
0175
0176
0177 void DicomFileMgr::AddMaterial(std::vector<G4String> wl)
0178 {
0179 if (theMaterials.size() > 0 && bMaterialsDensity) {
0180 G4Exception(
0181 "DicomFileMgr::AddMaterial", "DFM005", FatalException,
0182 "Trying to add a Material with :MATE and another with :MATE_DENS, check your input file");
0183 }
0184 bMaterialsDensity = false;
0185 theMaterials[G4UIcommand::ConvertToDouble(wl[2])] = wl[1];
0186 }
0187
0188
0189 void DicomFileMgr::AddMaterialDensity(std::vector<G4String> wl)
0190 {
0191 if (theMaterialsDensity.size() > 0 && !bMaterialsDensity) {
0192 G4Exception(
0193 "DicomFileMgr::AddMaterial", "DFM005", FatalException,
0194 "Trying to add a Material with :MATE and another with :MATE_DENS, check your input file");
0195 }
0196 bMaterialsDensity = true;
0197 theMaterialsDensity[G4UIcommand::ConvertToDouble(wl[2])] = wl[1];
0198 }
0199
0200
0201 void DicomFileMgr::AddCT2Density(std::vector<G4String> wl)
0202 {
0203 theCT2Density[G4UIcommand::ConvertToInt(wl[1])] = G4UIcommand::ConvertToDouble(wl[2]);
0204 G4cout << this << " AddCT2density " << theCT2Density.size() << G4endl;
0205 }
0206
0207
0208 G4double DicomFileMgr::Hounsfield2density(Uint32 Hval)
0209 {
0210 if (theCT2Density.size() == 0) {
0211 G4Exception("Hounsfield2density", "DCM004", FatalException, "No :CT2D line in input file");
0212 }
0213 std::map<G4int, G4double>::const_iterator ite = theCT2Density.begin();
0214 G4int minHval = (*ite).first;
0215 if (G4int(Hval) < minHval) {
0216 G4Exception("DicomHandler::Hounsfield2density", "", FatalException,
0217 ("Hval value too small, change input file " + std::to_string(Hval) + " < "
0218 + std::to_string(minHval))
0219 .c_str());
0220 }
0221
0222 ite = theCT2Density.end();
0223 ite--;
0224 G4int maxHval = (*ite).first;
0225 if (G4int(Hval) > maxHval) {
0226 G4Exception("DicomHandler::Hval2density", "", FatalException,
0227 ("Hval value too big, change CT2Density.dat file " + std::to_string(Hval) + " > "
0228 + std::to_string(maxHval))
0229 .c_str());
0230 }
0231
0232 G4float density = -1.;
0233 G4double deltaCT = 0;
0234 G4double deltaDensity = 0;
0235
0236 ite = theCT2Density.upper_bound(Hval);
0237 std::map<G4int, G4double>::const_iterator itePrev = ite;
0238 itePrev--;
0239
0240 deltaCT = (*ite).first - (*itePrev).first;
0241 deltaDensity = (*ite).second - (*itePrev).second;
0242
0243
0244 density = (*ite).second - (((*ite).first - Hval) * deltaDensity / deltaCT);
0245
0246 if (density < 0.) {
0247 G4Exception("DicomFileMgr::Hounsfiled2Density", "DFM002", FatalException,
0248 G4String("@@@ Error negative density = " + std::to_string(density)
0249 + " from HV = " + std::to_string(Hval))
0250 .c_str());
0251 }
0252
0253
0254 return density;
0255 }
0256
0257
0258 size_t DicomFileMgr::GetMaterialIndex(G4double Hval)
0259 {
0260 std::map<G4double, G4String>::iterator ite = theMaterials.upper_bound(Hval);
0261 if (ite == theMaterials.end()) {
0262 ite--;
0263 G4Exception("DicomFileMgr::GetMaterialIndex", "DFM004", FatalException,
0264 ("Hounsfiled value too big, change input file " + std::to_string(Hval) + " > "
0265 + std::to_string((*ite).first))
0266 .c_str());
0267 }
0268
0269 size_t dist = std::distance(theMaterials.begin(), ite);
0270
0271 return dist;
0272 }
0273
0274
0275 size_t DicomFileMgr::GetMaterialIndexByDensity(G4double density)
0276 {
0277 std::map<G4double, G4String>::iterator ite = theMaterialsDensity.upper_bound(density);
0278 if (ite == theMaterialsDensity.end()) {
0279 ite--;
0280 G4Exception("DicomFileMgr::GetMaterialIndexByDensity", "DFM003", FatalException,
0281 ("Density too big, change input file " + std::to_string(density) + " > "
0282 + std::to_string((*ite).first))
0283 .c_str());
0284 }
0285
0286 size_t dist = std::distance(theMaterialsDensity.begin(), ite);
0287
0288 return dist;
0289 }
0290
0291
0292 void DicomFileMgr::ProcessFiles()
0293 {
0294 if (theCTFiles.size() == 0) {
0295 G4Exception("CheckCTSlices", "DCM004", JustWarning, "No :FILE of type CT in input file");
0296 }
0297 else {
0298 CheckCTSlices();
0299
0300 BuildCTMaterials();
0301
0302 MergeCTFiles();
0303 }
0304
0305 G4cout << " PROCESSING PET FILES " << thePETFiles.size() << G4endl;
0306 if (thePETFiles.size() != 0) {
0307 CheckPETSlices();
0308
0309 BuildPETActivities();
0310
0311 MergePETFiles();
0312 }
0313
0314 DumpToTextFile();
0315 }
0316
0317
0318 void DicomFileMgr::CheckCTSlices()
0319 {
0320 size_t nSlices = theCTFiles.size();
0321 G4cout << " DicomFileMgr::Checking CT slices: " << nSlices << G4endl;
0322
0323 G4bool uniformSliceThickness = true;
0324
0325 if (nSlices > 1) {
0326 if (nSlices == 2) {
0327 mdct::const_iterator ite = theCTFiles.begin();
0328 DicomFileCT* one = (*ite).second;
0329 ite++;
0330 DicomFileCT* two = (*ite).second;
0331
0332 G4double real_distance = (two->GetLocation() - one->GetLocation()) / 2.;
0333
0334 if (one->GetMaxZ() != two->GetMinZ()) {
0335 one->SetMaxZ(one->GetLocation() + real_distance);
0336 two->SetMinZ(two->GetLocation() - real_distance);
0337
0338
0339 if (uniformSliceThickness) {
0340 one->SetMinZ(one->GetLocation() - real_distance);
0341 two->SetMaxZ(two->GetLocation() + real_distance);
0342 }
0343 }
0344 }
0345 else {
0346 mdct::iterator ite0 = theCTFiles.begin();
0347 mdct::iterator ite1 = ite0;
0348 ite1++;
0349 mdct::iterator ite2 = ite1;
0350 ite2++;
0351 for (; ite2 != theCTFiles.end(); ++ite0, ++ite1, ++ite2) {
0352 DicomFileCT* prev = (DicomFileCT*)((*ite0).second);
0353 DicomFileCT* slice = (DicomFileCT*)((*ite1).second);
0354 DicomFileCT* next = (DicomFileCT*)((*ite2).second);
0355 G4double real_up_distance = (next->GetLocation() - slice->GetLocation()) / 2.;
0356 G4double real_down_distance = (slice->GetLocation() - prev->GetLocation()) / 2.;
0357 G4double real_distance = real_up_distance + real_down_distance;
0358 G4double stated_distance = slice->GetMaxZ() - slice->GetMinZ();
0359
0360 if (std::fabs(real_distance - stated_distance) > 1.E-9) {
0361 unsigned int sliceNum = std::distance(theCTFiles.begin(), ite1);
0362 G4cerr << "\tDicomFileMgr::CheckCTSlices - Slice Distance Error in slice [" << sliceNum
0363 << "]: Distance between this slice and slices up and down = " << real_distance
0364 << " <> Slice width = " << stated_distance << " Slice locations "
0365 << prev->GetLocation() << " : " << slice->GetLocation() << " : "
0366 << next->GetLocation() << " DIFFERENCE= " << real_distance - stated_distance
0367 << G4endl;
0368 G4cerr << "!! WARNING: Geant4 will reset slice width so that it extends between "
0369 << "lower and upper slice " << G4endl;
0370
0371 slice->SetMinZ(slice->GetLocation() - real_down_distance);
0372 slice->SetMaxZ(slice->GetLocation() + real_up_distance);
0373
0374 if (ite0 == theCTFiles.begin()) {
0375 prev->SetMaxZ(slice->GetMinZ());
0376
0377
0378 if (uniformSliceThickness) {
0379 prev->SetMinZ(prev->GetLocation() - real_down_distance);
0380 }
0381 }
0382 if (static_cast<unsigned int>(std::distance(theCTFiles.begin(), ite2) + 1) == nSlices) {
0383 next->SetMinZ(slice->GetMaxZ());
0384
0385
0386 if (uniformSliceThickness) {
0387 next->SetMaxZ(next->GetLocation() + real_up_distance);
0388 }
0389 }
0390 }
0391 }
0392 }
0393 }
0394 }
0395
0396
0397 void DicomFileMgr::CheckPETSlices()
0398 {
0399 size_t nSlices = thePETFiles.size();
0400 G4cout << " DicomFileMgr::Checking PET slices: " << nSlices << G4endl;
0401
0402 G4bool uniformSliceThickness = true;
0403
0404 if (nSlices > 1) {
0405 if (nSlices == 2) {
0406 mdpet::const_iterator ite = thePETFiles.begin();
0407 DicomFilePET* one = (*ite).second;
0408 ite++;
0409 DicomFilePET* two = (*ite).second;
0410
0411 G4double real_distance = (two->GetLocation() - one->GetLocation()) / 2.;
0412
0413 if (one->GetMaxZ() != two->GetMinZ()) {
0414 one->SetMaxZ(one->GetLocation() + real_distance);
0415 two->SetMinZ(two->GetLocation() - real_distance);
0416
0417
0418 if (uniformSliceThickness) {
0419 one->SetMinZ(one->GetLocation() - real_distance);
0420 two->SetMaxZ(two->GetLocation() + real_distance);
0421 }
0422 }
0423 }
0424 else {
0425 mdpet::iterator ite0 = thePETFiles.begin();
0426 mdpet::iterator ite1 = ite0;
0427 ite1++;
0428 mdpet::iterator ite2 = ite1;
0429 ite2++;
0430 for (; ite2 != thePETFiles.end(); ++ite0, ++ite1, ++ite2) {
0431 DicomFilePET* prev = (DicomFilePET*)((*ite0).second);
0432 DicomFilePET* slice = (DicomFilePET*)((*ite1).second);
0433 DicomFilePET* next = (DicomFilePET*)((*ite2).second);
0434 G4double real_up_distance = (next->GetLocation() - slice->GetLocation()) / 2.;
0435 G4double real_down_distance = (slice->GetLocation() - prev->GetLocation()) / 2.;
0436 G4double real_distance = real_up_distance + real_down_distance;
0437 G4double stated_distance = slice->GetMaxZ() - slice->GetMinZ();
0438
0439 if (std::fabs(real_distance - stated_distance) > 1.E-9) {
0440 unsigned int sliceNum = std::distance(thePETFiles.begin(), ite1);
0441 G4cerr << "\tDicomFileMgr::CheckPETSlices - Slice Distance Error in slice [" << sliceNum
0442 << "]: Distance between this slice and slices up and down = " << real_distance
0443 << " <> Slice width = " << stated_distance << " Slice locations "
0444 << prev->GetLocation() << " : " << slice->GetLocation() << " : "
0445 << next->GetLocation() << " DIFFERENCE= " << real_distance - stated_distance
0446 << G4endl;
0447 G4cerr << "!! WARNING: Geant4 will reset slice width so that it extends between "
0448 << "lower and upper slice " << G4endl;
0449
0450 slice->SetMinZ(slice->GetLocation() - real_down_distance);
0451 slice->SetMaxZ(slice->GetLocation() + real_up_distance);
0452
0453 if (ite0 == thePETFiles.begin()) {
0454 prev->SetMaxZ(slice->GetMinZ());
0455
0456
0457 if (uniformSliceThickness) {
0458 prev->SetMinZ(prev->GetLocation() - real_down_distance);
0459 }
0460 }
0461 if (static_cast<unsigned int>(std::distance(thePETFiles.begin(), ite2) + 1) == nSlices) {
0462 next->SetMinZ(slice->GetMaxZ());
0463
0464
0465 if (uniformSliceThickness) {
0466 next->SetMaxZ(next->GetLocation() + real_up_distance);
0467 }
0468 }
0469 }
0470 }
0471 }
0472 }
0473 }
0474
0475
0476 void DicomFileMgr::BuildCTMaterials()
0477 {
0478 G4cout << " DicomFileMgr::Building Materials " << theCTFiles.size() << G4endl;
0479 mdct::const_iterator ite = theCTFiles.begin();
0480 for (; ite != theCTFiles.end(); ite++) {
0481 (*ite).second->BuildMaterials();
0482 }
0483 }
0484
0485
0486 void DicomFileMgr::BuildPETActivities()
0487 {
0488 G4cout << " DicomFileMgr::Building PETData " << thePETFiles.size() << G4endl;
0489 mdpet::const_iterator ite = thePETFiles.begin();
0490 for (; ite != thePETFiles.end(); ite++) {
0491 (*ite).second->BuildActivities();
0492 }
0493 }
0494
0495
0496 void DicomFileMgr::MergeCTFiles()
0497 {
0498 G4cout << " DicomFileMgr::Merging CT Files " << theCTFiles.size() << G4endl;
0499 mdct::const_iterator ite = theCTFiles.begin();
0500 theCTFileAll = new DicomFileCT(*((*ite).second));
0501 ite++;
0502 for (; ite != theCTFiles.end(); ite++) {
0503 (*theCTFileAll) += *((*ite).second);
0504 }
0505 }
0506
0507
0508 void DicomFileMgr::MergePETFiles()
0509 {
0510 G4cout << " DicomFileMgr::Merging PET Files " << thePETFiles.size() << G4endl;
0511 mdpet::const_iterator ite = thePETFiles.begin();
0512 thePETFileAll = new DicomFilePET(*((*ite).second));
0513 ite++;
0514 for (; ite != thePETFiles.end(); ite++) {
0515 (*thePETFileAll) += *((*ite).second);
0516 }
0517 }
0518
0519
0520 void DicomFileMgr::DumpToTextFile()
0521 {
0522 G4cout << " DicomFileMgr::Dumping To Text File " << G4endl;
0523 if (theCTFiles.size() != 0) {
0524 std::ofstream fout(theFileOutName);
0525
0526 if (!bMaterialsDensity) {
0527 fout << theMaterials.size() << std::endl;
0528 std::map<G4double, G4String>::const_iterator ite;
0529 G4int ii = 0;
0530 for (ite = theMaterials.begin(); ite != theMaterials.end(); ite++, ii++) {
0531 fout << ii << " \"" << (*ite).second << "\"" << std::endl;
0532 }
0533 }
0534 else {
0535 fout << theMaterialsDensity.size() << std::endl;
0536 std::map<G4double, G4String>::const_iterator ite;
0537 G4int ii = 0;
0538 for (ite = theMaterialsDensity.begin(); ite != theMaterialsDensity.end(); ite++, ii++) {
0539 fout << ii << " \"" << (*ite).second << "\"" << std::endl;
0540 }
0541 }
0542
0543 theCTFileAll->DumpHeaderToTextFile(fout);
0544 for (mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++) {
0545 (*itect).second->DumpMateIDsToTextFile(fout);
0546 }
0547 for (mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++) {
0548 (*itect).second->DumpDensitiesToTextFile(fout);
0549 }
0550 for (mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++) {
0551 (*itect).second->BuildStructureIDs();
0552 (*itect).second->DumpStructureIDsToTextFile(fout);
0553 }
0554
0555 std::vector<DicomFileStructure*> dfs = GetStructFiles();
0556 for (size_t i1 = 0; i1 < dfs.size(); i1++) {
0557 std::vector<DicomROI*> rois = dfs[i1]->GetROIs();
0558 for (size_t i2 = 0; i2 < rois.size(); i2++) {
0559 fout << rois[i2]->GetNumber() + 1 << " \"" << rois[i2]->GetName() << "\"" << G4endl;
0560 }
0561 }
0562 }
0563
0564 if (thePETFiles.size() != 0) {
0565 std::ofstream fout(theFileOutName);
0566
0567 thePETFileAll->DumpHeaderToTextFile(fout);
0568 for (mdpet::const_iterator itect = thePETFiles.begin(); itect != thePETFiles.end(); itect++) {
0569 (*itect).second->DumpActivitiesToTextFile(fout);
0570 }
0571 }
0572
0573 for (size_t i1 = 0; i1 < thePlanFiles.size(); i1++) {
0574 thePlanFiles[i1]->DumpToFile();
0575 }
0576 }
