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 // 27 /// \file medical/DICOM/src/DicomPartialDetect 28 /// \brief Implementation of the DicomPartialD 29 // 30 31 #include "DicomPartialDetectorConstruction.hh" 32 33 #include "G4Box.hh" 34 #include "G4Colour.hh" 35 #include "G4Element.hh" 36 #include "G4LogicalVolume.hh" 37 #include "G4Material.hh" 38 #include "G4NistManager.hh" 39 #include "G4PVParameterised.hh" 40 #include "G4PVPlacement.hh" 41 #include "G4PartialPhantomParameterisation.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4Tubs.hh" 44 #include "G4UIcommand.hh" 45 #include "G4VPhysicalVolume.hh" 46 #include "G4VisAttributes.hh" 47 #include "G4ios.hh" 48 #include "G4tgbMaterialMgr.hh" 49 #include "G4tgbVolumeMgr.hh" 50 #include "G4tgrMessenger.hh" 51 #include "globals.hh" 52 53 //....oooOO0OOooo........oooOO0OOooo........oo 54 DicomPartialDetectorConstruction::DicomPartial 55 : DicomDetectorConstruction(), 56 fPartialPhantomParam(0), 57 fNoVoxels(0), 58 fDimX(0), 59 fDimY(0), 60 fDimZ(0), 61 fOffsetX(0), 62 fOffsetY(0), 63 fOffsetZ(0) 64 {} 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 DicomPartialDetectorConstruction::~DicomPartia 68 69 //....oooOO0OOooo........oooOO0OOooo........oo 70 G4VPhysicalVolume* DicomPartialDetectorConstru 71 { 72 InitialisationOfMaterials(); 73 74 //----- Build world 75 G4double worldXDimension = 1. * m; 76 G4double worldYDimension = 1. * m; 77 G4double worldZDimension = 1. * m; 78 79 G4Box* world_box = new G4Box("WorldSolid", w 80 81 fWorld_logic = new G4LogicalVolume(world_box 82 83 G4VPhysicalVolume* world_phys = 84 new G4PVPlacement(0, G4ThreeVector(0, 0, 0 85 86 ReadPhantomData(); 87 88 ConstructPhantomContainer(); 89 ConstructPhantom(); 90 91 return world_phys; 92 } 93 94 //....oooOO0OOooo........oooOO0OOooo........oo 95 void DicomPartialDetectorConstruction::ReadPha 96 { 97 G4String dataFile = "Data.dat"; 98 std::ifstream finDF(dataFile.c_str()); 99 G4String fname; 100 if (finDF.good() != 1) { 101 G4Exception(" DicomPartialDetectorConstruc 102 " Problem reading data file: D 103 } 104 105 G4int compression; 106 finDF >> compression; // not used here 107 108 G4int numFiles; 109 finDF >> numFiles; // only 1 file supported 110 if (numFiles != 1) { 111 G4Exception("DicomPartialDetectorConstruct 112 "More than 1 DICOM file is not 113 } 114 for (G4int i = 0; i < numFiles; i++) { 115 finDF >> fname; 116 //--- Read one data file 117 fname += ".g4pdcm"; 118 ReadPhantomDataFile(fname); 119 } 120 121 finDF.close(); 122 } 123 124 //....oooOO0OOooo........oooOO0OOooo........oo 125 void DicomPartialDetectorConstruction::ReadPha 126 { 127 // G4String filename = "phantom.g4pdcm"; 128 #ifdef G4VERBOSE 129 G4cout << " DicomDetectorConstruction::ReadP 130 #endif 131 std::ifstream fin(fname.c_str(), std::ios_ba 132 if (!fin.is_open()) { 133 G4Exception("DicomPartialDetectorConstruct 134 G4String("File not found " + f 135 } 136 G4int nMaterials; 137 fin >> nMaterials; 138 G4String stemp; 139 G4int nmate; 140 for (G4int ii = 0; ii < nMaterials; ii++) { 141 fin >> nmate >> stemp; 142 G4cout << "DicomPartialDetectorConstructio 143 << nmate << " mate " << stemp << G4 144 if (ii != nmate) 145 G4Exception("DicomPartialDetectorConstru 146 "Material number should be i 147 wrong material numbe 148 } 149 150 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxels 151 G4cout << "DicomPartialDetectorConstruction: 152 << fNoVoxelsY << " " << fNoVoxelsZ << 153 fin >> fOffsetX >> fDimX; 154 fDimX = (fDimX - fOffsetX) / fNoVoxelsX; 155 fin >> fOffsetY >> fDimY; 156 fDimY = (fDimY - fOffsetY) / fNoVoxelsY; 157 fin >> fOffsetZ >> fDimZ; 158 fDimZ = (fDimZ - fOffsetZ) / fNoVoxelsZ; 159 G4cout << "DicomPartialDetectorConstruction: 160 << fOffsetX << G4endl; 161 G4cout << "DicomPartialDetectorConstruction: 162 << fOffsetY << G4endl; 163 G4cout << "DicomPartialDetectorConstruction: 164 << fOffsetZ << G4endl; 165 166 //--- Read voxels that are filled 167 fNoVoxels = 0; 168 // G4bool* isFilled = new G4bool[fNoVoxelsX 169 // fFilledIDs = new size_t[fNoVoxelsZ*fNoVo 170 //? fFilledIDs.insert(0); 171 G4int ifxmin1, ifxmax1; 172 for (G4int iz = 0; iz < fNoVoxelsZ; iz++) { 173 std::map<G4int, G4int> ifmin, ifmax; 174 for (G4int iy = 0; iy < fNoVoxelsY; iy++) 175 fin >> ifxmin1 >> ifxmax1; 176 // check coherence ... 177 178 ifmin[iy] = ifxmin1; 179 ifmax[iy] = ifxmax1; 180 G4int ifxdiff = ifxmax1 - ifxmin1 + 1; 181 if (ifxmax1 == -1 && ifxmin1 == -1) ifxd 182 fFilledIDs.insert(std::pair<G4int, G4int 183 G4cout << "DicomPartialDetectorConstruct 184 << " FilledIDs " << ifxdiff + fNo 185 << " N= " << fFilledIDs.size() << 186 // filledIDs[iz*fNoVoxelsY+iy+1] = ifxma 187 for (G4int ix = 0; ix < fNoVoxelsX; ix++ 188 if (ix >= G4int(ifxmin1) && ix <= G4in 189 fNoVoxels++; 190 } 191 else { 192 } 193 } 194 } 195 fFilledMins[iz] = ifmin; 196 fFilledMaxs[iz] = ifmax; 197 } 198 199 //--- Read material IDs 200 G4int mateID1; 201 fMateIDs = new size_t[fNoVoxelsX * fNoVoxels 202 G4int copyNo = 0; 203 for (G4int iz = 0; iz < fNoVoxelsZ; iz++) { 204 std::map<G4int, G4int> ifmin = fFilledMins 205 std::map<G4int, G4int> ifmax = fFilledMaxs 206 for (G4int iy = 0; iy < fNoVoxelsY; iy++) 207 ifxmin1 = ifmin[iy]; 208 ifxmax1 = ifmax[iy]; 209 for (G4int ix = 0; ix < fNoVoxelsX; ix++ 210 if (ix >= G4int(ifxmin1) && ix <= G4in 211 fin >> mateID1; 212 if (mateID1 < 0 || mateID1 >= nMater 213 G4Exception("DicomPartialDetectorC 214 G4String("Wrong index 215 + G4UIcommand 216 + G4UIcommand 217 .c_str()); 218 } 219 fMateIDs[copyNo] = mateID1; 220 copyNo++; 221 } 222 } 223 } 224 } 225 226 ReadVoxelDensitiesPartial(fin); 227 228 fin.close(); 229 } 230 231 //....oooOO0OOooo........oooOO0OOooo........oo 232 void DicomPartialDetectorConstruction::ReadVox 233 { 234 std::map<G4int, std::pair<G4double, G4double 235 std::map<G4int, std::pair<G4double, G4double 236 for (G4int ii = 0; ii < G4int(fOriginalMater 237 densiMinMax[ii] = std::pair<G4double, G4do 238 } 239 240 //----- Define density differences (maximum 241 // a new material) 242 char* part = std::getenv("DICOM_CHANGE_MATER 243 G4double densityDiff = -1.; 244 if (part) densityDiff = G4UIcommand::Convert 245 std::map<G4int, G4double> densityDiffs; 246 if (densityDiff != -1.) { 247 for (G4int ii = 0; ii < G4int(fOriginalMat 248 densityDiffs[ii] = densityDiff; // curr 249 // with same step 250 } 251 } 252 else { 253 if (fMaterials.size() == 0) { // do it on 254 for (unsigned int ii = 0; ii < fOriginal 255 fMaterials.push_back(fOriginalMaterial 256 } 257 } 258 } 259 // densityDiffs[0] = 0.0001; //fAir 260 261 //--- Calculate the average material density 262 std::map<std::pair<G4Material*, G4int>, matI 263 264 G4double dens1; 265 G4int ifxmin1, ifxmax1; 266 267 //---- Read the material densities 268 G4int copyNo = 0; 269 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) { 270 std::map<G4int, G4int> ifmin = fFilledMins 271 std::map<G4int, G4int> ifmax = fFilledMaxs 272 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) 273 ifxmin1 = ifmin[iy]; 274 ifxmax1 = ifmax[iy]; 275 for (G4int ix = 0; ix < fNoVoxelsX; ++ix 276 if (ix >= G4int(ifxmin1) && ix <= G4in 277 fin >> dens1; 278 //--- store the minimum and maximum 279 //(just for printing) 280 mpite = densiMinMax.find(G4int(fMate 281 if (dens1 < (*mpite).second.first) ( 282 if (dens1 > (*mpite).second.second) 283 284 //--- Get material from original lis 285 G4int mateID = G4int(fMateIDs[copyNo 286 G4Material* mate = fOriginalMaterial 287 // G4cout << copyNo << " mate 288 //--- Check if density is equal to t 289 if (std::fabs(dens1 - mate->GetDensi 290 291 //--- Build material name with fOrig 292 // float densityBin = density 293 // (G4int(dens1/densityDiffs[ 294 G4String newMateName = mate->GetName 295 296 G4int densityBin = 0; 297 // G4cout << " densityBin " < 298 // << dens1 << " " <<densityD 299 300 if (densityDiff != -1.) { 301 densityBin = (G4int(dens1 / densit 302 newMateName = mate->GetName() + G4 303 } 304 305 //--- Look if it is the first voxel 306 std::pair<G4Material*, G4int> matden 307 308 auto mppite = newMateDens.find(matde 309 if (mppite != newMateDens.cend()) { 310 matInfo* mi = (*mppite).second; 311 mi->fSumdens += dens1; 312 mi->fNvoxels++; 313 fMateIDs[copyNo] = fOriginalMateri 314 // G4cout << copyNo << " mat new 315 //<< fOriginalMaterials.size()-1 + 316 } 317 else { 318 matInfo* mi = new matInfo; 319 mi->fSumdens = dens1; 320 mi->fNvoxels = 1; 321 mi->fId = G4int(newMateDens.size() 322 newMateDens[matdens] = mi; 323 fMateIDs[copyNo] = fOriginalMateri 324 // G4cout << copyNo << " 325 // << fOriginalMaterials. 326 } 327 copyNo++; 328 // G4cout << ix << " " << iy 329 // << " filling fMateIDs " << copyN 330 // fMateIDs[copyNo] = atoi(ci 331 } 332 } 333 } 334 } 335 336 //----- Build the list of phantom materials 337 //--- Add original materials 338 for (auto mimite = fOriginalMaterials.cbegin 339 fPhantomMaterials.push_back((*mimite)); 340 } 341 // 342 //---- Build and add new materials 343 for (auto mppite = newMateDens.cbegin(); mpp 344 G4double averdens = (*mppite).second->fSum 345 G4double saverdens = G4int(1000.001 * aver 346 G4cout << "GmReadPhantomGeometry::ReadVoxe 347 << saverdens << " -> " << G4int(100 348 << " " << G4int(1000 * averdens) / 349 350 G4cout << "GmReadPhantomGeometry::ReadVoxe 351 << saverdens << " -> " << G4UIcomma 352 << (*mppite).second->fNvoxels << G4 353 G4String mateName = 354 ((*mppite).first).first->GetName() + "_" 355 fPhantomMaterials.push_back( 356 BuildMaterialWithChangingDensity((*mppit 357 } 358 } 359 360 //....oooOO0OOooo........oooOO0OOooo........oo 361 void DicomPartialDetectorConstruction::Constru 362 { 363 // build a clinder that encompass the partia 364 // /dicom/intersectWithUserVolume 0. 0. 0. 365 //---- Extract number of voxels and voxel di 366 367 G4String contname = "PhantomContainer"; 368 369 //----- Define the volume that contains all 370 G4Tubs* container_tube = new G4Tubs(contname 371 372 fContainer_logic = new G4LogicalVolume(conta 373 374 G4ThreeVector posCentreVoxels(0., 0., 0.); 375 // G4cout << " PhantomContainer posCentreVox 376 G4RotationMatrix* rotm = new G4RotationMatri 377 378 fContainer_phys = new G4PVPlacement(rotm, / 379 posCentr 380 fContain 381 contname 382 fWorld_l 383 false, 384 1); // 385 } 386 387 //....oooOO0OOooo........oooOO0OOooo........oo 388 void DicomPartialDetectorConstruction::Constru 389 { 390 G4String OptimAxis = "kXAxis"; 391 G4bool bSkipEqualMaterials = 0; 392 G4int regStructureID = 2; 393 394 G4ThreeVector posCentreVoxels(0., 0., 0.); 395 396 //----- Build phantom 397 G4String voxelName = "phantom"; 398 G4Box* phantom_solid = new G4Box(voxelName, 399 G4LogicalVolume* phantom_logic = new G4Logic 400 G4bool pVis = 0; 401 if (!pVis) { 402 G4VisAttributes* visAtt = new G4VisAttribu 403 phantom_logic->SetVisAttributes(visAtt); 404 } 405 406 G4double theSmartless = 0.5; 407 // fContainer_logic->SetSmartless( theSmart 408 phantom_logic->SetSmartless(theSmartless); 409 410 fPartialPhantomParam = new G4PartialPhantomP 411 412 fPartialPhantomParam->SetMaterials(fPhantomM 413 fPartialPhantomParam->SetVoxelDimensions(fDi 414 fPartialPhantomParam->SetNoVoxels(fNoVoxelsX 415 fPartialPhantomParam->SetMaterialIndices(fMa 416 417 fPartialPhantomParam->SetFilledIDs(fFilledID 418 419 fPartialPhantomParam->SetFilledMins(fFilledM 420 421 fPartialPhantomParam->BuildContainerWalls(); 422 423 // G4cout << " Number of Materials " << fPh 424 // G4cout << " SetMaterialIndices(0) " << f 425 426 G4PVParameterised* phantom_phys = 0; 427 if (OptimAxis == "kUndefined") { 428 phantom_phys = new G4PVParameterised(voxel 429 fNoVo 430 } 431 else if (OptimAxis == "kXAxis") { 432 // G4cout << " optim kX " << G4endl; 433 phantom_phys = new G4PVParameterised(voxel 434 fNoVo 435 fPartialPhantomParam->SetSkipEqualMaterial 436 } 437 else { 438 G4Exception("GmReadPhantomGeometry::Constr 439 G4String("Wrong argument to pa 440 GmReadPhantomGeometry:Phantom:OptimAx 441 Only allowed 'kUndefined' or 'kXAxis' 442 + OptimAxis) 443 .c_str()); 444 } 445 446 phantom_phys->SetRegularStructureId(regStruc 447 } 448