Warning, file /geant4/examples/extended/medical/DICOM/G4DicomReader/src/DicomFileMgr.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
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 }
