Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 #include "DicomFileMgr.hh" 27 28 #include "DicomFileCT.hh" 29 #include "DicomFilePET.hh" 30 #include "DicomFilePlan.hh" 31 #include "DicomFileStructure.hh" 32 #include "dcmtk/dcmdata/dcdeftag.h" 33 34 #include "G4UIcommand.hh" 35 #include "G4tgrFileIn.hh" 36 37 DicomFileMgr* DicomFileMgr::theInstance = 0; 38 int DicomFileMgr::verbose = 1; 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 DicomFileMgr* DicomFileMgr::GetInstance() 42 { 43 if (!theInstance) { 44 theInstance = new DicomFileMgr; 45 } 46 return theInstance; 47 } 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 DicomFileMgr::DicomFileMgr() 51 { 52 fCompression = 1.; 53 theCTFileAll = 0; 54 theStructureNCheck = 4; 55 theStructureNMaxROI = 100; 56 } 57 58 //....oooOO0OOooo........oooOO0OOooo........oo 59 void DicomFileMgr::Convert(G4String fileName) 60 { 61 G4tgrFileIn fin = G4tgrFileIn::GetInstance(f 62 std::vector<G4String> wl; 63 // Read each file in file list 64 theFileOutName = "test.g4dcm"; 65 int ii; 66 for (ii = 0;; ii++) { 67 if (!fin.GetWordsInLine(wl)) break; 68 if (wl[0] == ":COMPRESSION") { 69 CheckNColumns(wl, 2); 70 SetCompression(wl[1]); 71 } 72 else if (wl[0] == ":FILE") { 73 CheckNColumns(wl, 2); 74 G4cout << "@@@@@@@ Reading FILE: " << wl 75 AddFile(wl[1]); 76 } 77 else if (wl[0] == ":FILE_OUT") { 78 CheckNColumns(wl, 2); 79 theFileOutName = wl[1]; 80 } 81 else if (wl[0] == ":MATE_DENS") { 82 CheckNColumns(wl, 3); 83 AddMaterialDensity(wl); 84 } 85 else if (wl[0] == ":MATE") { 86 CheckNColumns(wl, 3); 87 AddMaterial(wl); 88 } 89 else if (wl[0] == ":CT2D") { 90 CheckNColumns(wl, 3); 91 AddCT2Density(wl); 92 } 93 else { 94 G4Exception("DICOM2G4", "Wrong argument" 95 G4String("UNKNOWN TAG IN FIL 96 } 97 } 98 99 //@@@@@@ Process files 100 ProcessFiles(); 101 } 102 103 //....oooOO0OOooo........oooOO0OOooo........oo 104 void DicomFileMgr::CheckNColumns(std::vector<G 105 { 106 if (wl.size() != vsizeTh) { 107 G4cerr << " Reading line " << G4endl; 108 for (size_t ii = 0; ii < wl.size(); ii++) 109 G4cerr << wl[ii] << " "; 110 } 111 G4cerr << G4endl; 112 G4Exception("DICOM2G4", "D2G0010", FatalEr 113 ("Wrong number of columns in l 114 + std::to_string(vsizeTh)) 115 .c_str()); 116 } 117 } 118 119 //....oooOO0OOooo........oooOO0OOooo........oo 120 void DicomFileMgr::SetCompression(G4String fCo 121 { 122 fCompression = G4UIcommand::ConvertToDouble( 123 } 124 125 //....oooOO0OOooo........oooOO0OOooo........oo 126 void DicomFileMgr::AddFile(G4String fileName) 127 { 128 DcmFileFormat dfile; 129 if (!(dfile.loadFile(fileName.c_str())).good 130 G4Exception("DicomHandler::ReadFile", "dfi 131 ("Error reading file " + fileN 132 } 133 DcmDataset* dset = dfile.getDataset(); 134 135 OFString dModality; 136 if (!dset->findAndGetOFString(DCM_Modality, 137 G4Exception("DicomHandler::ReadData ", "", 138 } 139 140 if (dModality == "CT") { 141 DicomFileCT* df = new DicomFileCT(dset); 142 df->ReadData(); 143 df->SetFileName(fileName); 144 // reorder by location 145 theCTFiles[df->GetMaxZ()] = df; 146 } 147 else if (dModality == "RTSTRUCT") { 148 DicomFileStructure* df = new DicomFileStru 149 df->ReadData(); 150 df->SetFileName(fileName); 151 // theFiles.insert(msd::value_type(dMod 152 theStructFiles.push_back(df); 153 } 154 else if (dModality == "RTPLAN") { 155 DicomFilePlan* df = new DicomFilePlan(dset 156 df->ReadData(); 157 df->SetFileName(fileName); 158 // theFiles.insert(msd::value_type(dMod 159 thePlanFiles.push_back(df); 160 } 161 else if (dModality == "PT") { 162 DicomFilePET* df = new DicomFilePET(dset); 163 df->ReadData(); 164 df->SetFileName(fileName); 165 // theFiles.insert(msd::value_type(dMod 166 thePETFiles[df->GetMaxZ()] = df; 167 // thePETFiles.push_back(df); 168 } 169 else { 170 G4Exception( 171 "DicomFileMgr::AddFIle()", "DFM001", Fat 172 (G4String("File is not of type CT or RTS 173 } 174 } 175 176 //....oooOO0OOooo........oooOO0OOooo........oo 177 void DicomFileMgr::AddMaterial(std::vector<G4S 178 { 179 if (theMaterials.size() > 0 && bMaterialsDen 180 G4Exception( 181 "DicomFileMgr::AddMaterial", "DFM005", F 182 "Trying to add a Material with :MATE and 183 } 184 bMaterialsDensity = false; 185 theMaterials[G4UIcommand::ConvertToDouble(wl 186 } 187 188 //....oooOO0OOooo........oooOO0OOooo........oo 189 void DicomFileMgr::AddMaterialDensity(std::vec 190 { 191 if (theMaterialsDensity.size() > 0 && !bMate 192 G4Exception( 193 "DicomFileMgr::AddMaterial", "DFM005", F 194 "Trying to add a Material with :MATE and 195 } 196 bMaterialsDensity = true; 197 theMaterialsDensity[G4UIcommand::ConvertToDo 198 } 199 200 //....oooOO0OOooo........oooOO0OOooo........oo 201 void DicomFileMgr::AddCT2Density(std::vector<G 202 { 203 theCT2Density[G4UIcommand::ConvertToInt(wl[1 204 G4cout << this << " AddCT2density " << theCT 205 } 206 207 //....oooOO0OOooo........oooOO0OOooo........oo 208 G4double DicomFileMgr::Hounsfield2density(Uint 209 { 210 if (theCT2Density.size() == 0) { 211 G4Exception("Hounsfield2density", "DCM004" 212 } 213 std::map<G4int, G4double>::const_iterator it 214 G4int minHval = (*ite).first; 215 if (G4int(Hval) < minHval) { 216 G4Exception("DicomHandler::Hounsfield2dens 217 ("Hval value too small, change 218 + std::to_string(minHval)) 219 .c_str()); 220 } 221 222 ite = theCT2Density.end(); 223 ite--; 224 G4int maxHval = (*ite).first; 225 if (G4int(Hval) > maxHval) { 226 G4Exception("DicomHandler::Hval2density", 227 ("Hval value too big, change C 228 + std::to_string(maxHval)) 229 .c_str()); 230 } 231 232 G4float density = -1.; 233 G4double deltaCT = 0; 234 G4double deltaDensity = 0; 235 236 ite = theCT2Density.upper_bound(Hval); 237 std::map<G4int, G4double>::const_iterator it 238 itePrev--; 239 240 deltaCT = (*ite).first - (*itePrev).first; 241 deltaDensity = (*ite).second - (*itePrev).se 242 243 // interpolating linearly 244 density = (*ite).second - (((*ite).first - H 245 246 if (density < 0.) { 247 G4Exception("DicomFileMgr::Hounsfiled2Dens 248 G4String("@@@ Error negative d 249 + " from HV = " + std 250 .c_str()); 251 } 252 253 // G4cout << " Hval2density " << Hval << " 254 return density; 255 } 256 257 //....oooOO0OOooo........oooOO0OOooo........oo 258 size_t DicomFileMgr::GetMaterialIndex(G4double 259 { 260 std::map<G4double, G4String>::iterator ite = 261 if (ite == theMaterials.end()) { 262 ite--; 263 G4Exception("DicomFileMgr::GetMaterialInde 264 ("Hounsfiled value too big, ch 265 + std::to_string((*ite).first 266 .c_str()); 267 } 268 269 size_t dist = std::distance(theMaterials.beg 270 271 return dist; 272 } 273 274 //....oooOO0OOooo........oooOO0OOooo........oo 275 size_t DicomFileMgr::GetMaterialIndexByDensity 276 { 277 std::map<G4double, G4String>::iterator ite = 278 if (ite == theMaterialsDensity.end()) { 279 ite--; 280 G4Exception("DicomFileMgr::GetMaterialInde 281 ("Density too big, change inpu 282 + std::to_string((*ite).first 283 .c_str()); 284 } 285 286 size_t dist = std::distance(theMaterialsDens 287 288 return dist; 289 } 290 291 //....oooOO0OOooo........oooOO0OOooo........oo 292 void DicomFileMgr::ProcessFiles() 293 { 294 if (theCTFiles.size() == 0) { 295 G4Exception("CheckCTSlices", "DCM004", Jus 296 } 297 else { 298 CheckCTSlices(); 299 300 BuildCTMaterials(); 301 302 MergeCTFiles(); 303 } 304 305 G4cout << " PROCESSING PET FILES " << thePET 306 if (thePETFiles.size() != 0) { 307 CheckPETSlices(); 308 309 BuildPETActivities(); 310 311 MergePETFiles(); 312 } 313 314 DumpToTextFile(); 315 } 316 317 //....oooOO0OOooo........oooOO0OOooo........oo 318 void DicomFileMgr::CheckCTSlices() 319 { 320 size_t nSlices = theCTFiles.size(); 321 G4cout << " DicomFileMgr::Checking CT slices 322 323 G4bool uniformSliceThickness = true; 324 325 if (nSlices > 1) { 326 if (nSlices == 2) { 327 mdct::const_iterator ite = theCTFiles.be 328 DicomFileCT* one = (*ite).second; 329 ite++; 330 DicomFileCT* two = (*ite).second; 331 332 G4double real_distance = (two->GetLocati 333 334 if (one->GetMaxZ() != two->GetMinZ()) { 335 one->SetMaxZ(one->GetLocation() + real 336 two->SetMinZ(two->GetLocation() - real 337 // one->SetMinZ(one->GetLocation()-rea 338 // two->SetMaxZ(two->GetLocation()+rea 339 if (uniformSliceThickness) { 340 one->SetMinZ(one->GetLocation() - re 341 two->SetMaxZ(two->GetLocation() + re 342 } 343 } 344 } 345 else { 346 mdct::iterator ite0 = theCTFiles.begin() 347 mdct::iterator ite1 = ite0; 348 ite1++; 349 mdct::iterator ite2 = ite1; 350 ite2++; 351 for (; ite2 != theCTFiles.end(); ++ite0, 352 DicomFileCT* prev = (DicomFileCT*)((*i 353 DicomFileCT* slice = (DicomFileCT*)((* 354 DicomFileCT* next = (DicomFileCT*)((*i 355 G4double real_up_distance = (next->Get 356 G4double real_down_distance = (slice-> 357 G4double real_distance = real_up_dista 358 G4double stated_distance = slice->GetM 359 360 if (std::fabs(real_distance - stated_d 361 unsigned int sliceNum = std::distanc 362 G4cerr << "\tDicomFileMgr::CheckCTSl 363 << "]: Distance between this 364 << " <> Slice width = " << st 365 << prev->GetLocation() << " : 366 << next->GetLocation() << " D 367 << G4endl; 368 G4cerr << "!! WARNING: Geant4 will r 369 << "lower and upper slice " < 370 371 slice->SetMinZ(slice->GetLocation() 372 slice->SetMaxZ(slice->GetLocation() 373 374 if (ite0 == theCTFiles.begin()) { 375 prev->SetMaxZ(slice->GetMinZ()); 376 // Using below would make all slic 377 // prev->SetMinZ(prev->GetLocation 378 if (uniformSliceThickness) { 379 prev->SetMinZ(prev->GetLocation( 380 } 381 } 382 if (static_cast<unsigned int>(std::d 383 next->SetMinZ(slice->GetMaxZ()); 384 // Using below would make all slic 385 // next->SetMaxZ(next->GetLocation 386 if (uniformSliceThickness) { 387 next->SetMaxZ(next->GetLocation( 388 } 389 } 390 } 391 } 392 } 393 } 394 } 395 396 //....oooOO0OOooo........oooOO0OOooo........oo 397 void DicomFileMgr::CheckPETSlices() 398 { 399 size_t nSlices = thePETFiles.size(); 400 G4cout << " DicomFileMgr::Checking PET slice 401 402 G4bool uniformSliceThickness = true; 403 404 if (nSlices > 1) { 405 if (nSlices == 2) { 406 mdpet::const_iterator ite = thePETFiles. 407 DicomFilePET* one = (*ite).second; 408 ite++; 409 DicomFilePET* two = (*ite).second; 410 411 G4double real_distance = (two->GetLocati 412 413 if (one->GetMaxZ() != two->GetMinZ()) { 414 one->SetMaxZ(one->GetLocation() + real 415 two->SetMinZ(two->GetLocation() - real 416 // one->SetMinZ(one->GetLocation()-rea 417 // two->SetMaxZ(two->GetLocation()+rea 418 if (uniformSliceThickness) { 419 one->SetMinZ(one->GetLocation() - re 420 two->SetMaxZ(two->GetLocation() + re 421 } 422 } 423 } 424 else { 425 mdpet::iterator ite0 = thePETFiles.begin 426 mdpet::iterator ite1 = ite0; 427 ite1++; 428 mdpet::iterator ite2 = ite1; 429 ite2++; 430 for (; ite2 != thePETFiles.end(); ++ite0 431 DicomFilePET* prev = (DicomFilePET*)(( 432 DicomFilePET* slice = (DicomFilePET*)( 433 DicomFilePET* next = (DicomFilePET*)(( 434 G4double real_up_distance = (next->Get 435 G4double real_down_distance = (slice-> 436 G4double real_distance = real_up_dista 437 G4double stated_distance = slice->GetM 438 439 if (std::fabs(real_distance - stated_d 440 unsigned int sliceNum = std::distanc 441 G4cerr << "\tDicomFileMgr::CheckPETS 442 << "]: Distance between this 443 << " <> Slice width = " << st 444 << prev->GetLocation() << " : 445 << next->GetLocation() << " D 446 << G4endl; 447 G4cerr << "!! WARNING: Geant4 will r 448 << "lower and upper slice " < 449 450 slice->SetMinZ(slice->GetLocation() 451 slice->SetMaxZ(slice->GetLocation() 452 453 if (ite0 == thePETFiles.begin()) { 454 prev->SetMaxZ(slice->GetMinZ()); 455 // Using below would make all slic 456 // prev->SetMinZ(prev->GetLocation 457 if (uniformSliceThickness) { 458 prev->SetMinZ(prev->GetLocation( 459 } 460 } 461 if (static_cast<unsigned int>(std::d 462 next->SetMinZ(slice->GetMaxZ()); 463 // Using below would make all slic 464 // next->SetMaxZ(next->GetLocation 465 if (uniformSliceThickness) { 466 next->SetMaxZ(next->GetLocation( 467 } 468 } 469 } 470 } 471 } 472 } 473 } 474 475 //....oooOO0OOooo........oooOO0OOooo........oo 476 void DicomFileMgr::BuildCTMaterials() 477 { 478 G4cout << " DicomFileMgr::Building Materials 479 mdct::const_iterator ite = theCTFiles.begin( 480 for (; ite != theCTFiles.end(); ite++) { 481 (*ite).second->BuildMaterials(); 482 } 483 } 484 485 //....oooOO0OOooo........oooOO0OOooo........oo 486 void DicomFileMgr::BuildPETActivities() 487 { 488 G4cout << " DicomFileMgr::Building PETData " 489 mdpet::const_iterator ite = thePETFiles.begi 490 for (; ite != thePETFiles.end(); ite++) { 491 (*ite).second->BuildActivities(); 492 } 493 } 494 495 //....oooOO0OOooo........oooOO0OOooo........oo 496 void DicomFileMgr::MergeCTFiles() 497 { 498 G4cout << " DicomFileMgr::Merging CT Files " 499 mdct::const_iterator ite = theCTFiles.begin( 500 theCTFileAll = new DicomFileCT(*((*ite).seco 501 ite++; 502 for (; ite != theCTFiles.end(); ite++) { 503 (*theCTFileAll) += *((*ite).second); 504 } 505 } 506 507 //....oooOO0OOooo........oooOO0OOooo........oo 508 void DicomFileMgr::MergePETFiles() 509 { 510 G4cout << " DicomFileMgr::Merging PET Files 511 mdpet::const_iterator ite = thePETFiles.begi 512 thePETFileAll = new DicomFilePET(*((*ite).se 513 ite++; 514 for (; ite != thePETFiles.end(); ite++) { 515 (*thePETFileAll) += *((*ite).second); 516 } 517 } 518 519 //....oooOO0OOooo........oooOO0OOooo........oo 520 void DicomFileMgr::DumpToTextFile() 521 { 522 G4cout << " DicomFileMgr::Dumping To Text Fi 523 if (theCTFiles.size() != 0) { 524 std::ofstream fout(theFileOutName); 525 526 if (!bMaterialsDensity) { 527 fout << theMaterials.size() << std::endl 528 std::map<G4double, G4String>::const_iter 529 G4int ii = 0; 530 for (ite = theMaterials.begin(); ite != 531 fout << ii << " \"" << (*ite).second < 532 } 533 } 534 else { 535 fout << theMaterialsDensity.size() << st 536 std::map<G4double, G4String>::const_iter 537 G4int ii = 0; 538 for (ite = theMaterialsDensity.begin(); 539 fout << ii << " \"" << (*ite).second < 540 } 541 } 542 543 theCTFileAll->DumpHeaderToTextFile(fout); 544 for (mdct::const_iterator itect = theCTFil 545 (*itect).second->DumpMateIDsToTextFile(f 546 } 547 for (mdct::const_iterator itect = theCTFil 548 (*itect).second->DumpDensitiesToTextFile 549 } 550 for (mdct::const_iterator itect = theCTFil 551 (*itect).second->BuildStructureIDs(); 552 (*itect).second->DumpStructureIDsToTextF 553 } 554 555 std::vector<DicomFileStructure*> dfs = Get 556 for (size_t i1 = 0; i1 < dfs.size(); i1++) 557 std::vector<DicomROI*> rois = dfs[i1]->G 558 for (size_t i2 = 0; i2 < rois.size(); i2 559 fout << rois[i2]->GetNumber() + 1 << " 560 } 561 } 562 } 563 564 if (thePETFiles.size() != 0) { 565 std::ofstream fout(theFileOutName); 566 567 thePETFileAll->DumpHeaderToTextFile(fout); 568 for (mdpet::const_iterator itect = thePETF 569 (*itect).second->DumpActivitiesToTextFil 570 } 571 } 572 573 for (size_t i1 = 0; i1 < thePlanFiles.size() 574 thePlanFiles[i1]->DumpToFile(); 575 } 576 } 577