| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file medical/DICOM/src/DicomPartialDetectorConstruction.cc
28 /// \brief Implementation of the DicomPartialDetectorConstruction class
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........oooOO0OOooo........oooOO0OOooo......
54 DicomPartialDetectorConstruction::DicomPartialDetectorConstruction()
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........oooOO0OOooo........oooOO0OOooo......
67 DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction() {}
68
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 G4VPhysicalVolume* DicomPartialDetectorConstruction::Construct()
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", worldXDimension, worldYDimension, worldZDimension);
80
81 fWorld_logic = new G4LogicalVolume(world_box, fAir, "WorldLogical", 0, 0, 0);
82
83 G4VPhysicalVolume* world_phys =
84 new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "World", fWorld_logic, 0, false, 0);
85
86 ReadPhantomData();
87
88 ConstructPhantomContainer();
89 ConstructPhantom();
90
91 return world_phys;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 void DicomPartialDetectorConstruction::ReadPhantomData()
96 {
97 G4String dataFile = "Data.dat";
98 std::ifstream finDF(dataFile.c_str());
99 G4String fname;
100 if (finDF.good() != 1) {
101 G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
102 " Problem reading data file: Data.dat");
103 }
104
105 G4int compression;
106 finDF >> compression; // not used here
107
108 G4int numFiles;
109 finDF >> numFiles; // only 1 file supported for the moment
110 if (numFiles != 1) {
111 G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
112 "More than 1 DICOM file is not supported");
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........oooOO0OOooo........oooOO0OOooo......
125 void DicomPartialDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
126 {
127 // G4String filename = "phantom.g4pdcm";
128 #ifdef G4VERBOSE
129 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
130 #endif
131 std::ifstream fin(fname.c_str(), std::ios_base::in);
132 if (!fin.is_open()) {
133 G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile", "", FatalException,
134 G4String("File not found " + fname).c_str());
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 << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate " << ii << " = "
143 << nmate << " mate " << stemp << G4endl;
144 if (ii != nmate)
145 G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
146 "Material number should be in increasing order: \
147 wrong material number ");
148 }
149
150 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
151 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z " << fNoVoxelsX << " "
152 << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
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::ReadPhantomData voxelDimX " << fDimX << " fOffsetX "
160 << fOffsetX << G4endl;
161 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY " << fDimY << " fOffsetY "
162 << fOffsetY << G4endl;
163 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ " << fDimZ << " fOffsetZ "
164 << fOffsetZ << G4endl;
165
166 //--- Read voxels that are filled
167 fNoVoxels = 0;
168 // G4bool* isFilled = new G4bool[fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ];
169 // fFilledIDs = new size_t[fNoVoxelsZ*fNoVoxelsY+1];
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) ifxdiff = 0;
182 fFilledIDs.insert(std::pair<G4int, G4int>(ifxdiff + fNoVoxels - 1, ifxmin1));
183 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert "
184 << " FilledIDs " << ifxdiff + fNoVoxels - 1 << " min " << ifxmin1
185 << " N= " << fFilledIDs.size() << G4endl;
186 // filledIDs[iz*fNoVoxelsY+iy+1] = ifxmax1-ifxmin1+1;
187 for (G4int ix = 0; ix < fNoVoxelsX; ix++) {
188 if (ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1)) {
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 * fNoVoxelsY * fNoVoxelsZ];
202 G4int copyNo = 0;
203 for (G4int iz = 0; iz < fNoVoxelsZ; iz++) {
204 std::map<G4int, G4int> ifmin = fFilledMins[iz];
205 std::map<G4int, G4int> ifmax = fFilledMaxs[iz];
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 <= G4int(ifxmax1)) {
211 fin >> mateID1;
212 if (mateID1 < 0 || mateID1 >= nMaterials) {
213 G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalException,
214 G4String("Wrong index in phantom file: It should be between 0 and "
215 + G4UIcommand::ConvertToString(nMaterials - 1) + ", while it is "
216 + G4UIcommand::ConvertToString(mateID1))
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........oooOO0OOooo........oooOO0OOooo......
232 void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial(std::ifstream& fin)
233 {
234 std::map<G4int, std::pair<G4double, G4double>> densiMinMax;
235 std::map<G4int, std::pair<G4double, G4double>>::iterator mpite;
236 for (G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii) {
237 densiMinMax[ii] = std::pair<G4double, G4double>(DBL_MAX, -DBL_MAX);
238 }
239
240 //----- Define density differences (maximum density difference to create
241 // a new material)
242 char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
243 G4double densityDiff = -1.;
244 if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
245 std::map<G4int, G4double> densityDiffs;
246 if (densityDiff != -1.) {
247 for (G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii) {
248 densityDiffs[ii] = densityDiff; // currently all materials
249 // with same step
250 }
251 }
252 else {
253 if (fMaterials.size() == 0) { // do it only for first slice
254 for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
255 fMaterials.push_back(fOriginalMaterials[ii]);
256 }
257 }
258 }
259 // densityDiffs[0] = 0.0001; //fAir
260
261 //--- Calculate the average material density for each material/density bin
262 std::map<std::pair<G4Material*, G4int>, matInfo*> newMateDens;
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[iz];
271 std::map<G4int, G4int> ifmax = fFilledMaxs[iz];
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 <= G4int(ifxmax1)) {
277 fin >> dens1;
278 //--- store the minimum and maximum density for each material
279 //(just for printing)
280 mpite = densiMinMax.find(G4int(fMateIDs[copyNo]));
281 if (dens1 < (*mpite).second.first) (*mpite).second.first = dens1;
282 if (dens1 > (*mpite).second.second) (*mpite).second.second = dens1;
283
284 //--- Get material from original list of material in file
285 G4int mateID = G4int(fMateIDs[copyNo]);
286 G4Material* mate = fOriginalMaterials[mateID];
287 // G4cout << copyNo << " mateID " << mateID << G4endl;
288 //--- Check if density is equal to the original material density
289 if (std::fabs(dens1 - mate->GetDensity() / g * cm3) < 1.e-9) continue;
290
291 //--- Build material name with fOriginalMaterials name + density
292 // float densityBin = densityDiffs[mateID] *
293 // (G4int(dens1/densityDiffs[mateID])+0.5);
294 G4String newMateName = mate->GetName();
295
296 G4int densityBin = 0;
297 // G4cout << " densityBin " << densityBin << " "
298 // << dens1 << " " <<densityDiffs[mateID] << G4endl;
299
300 if (densityDiff != -1.) {
301 densityBin = (G4int(dens1 / densityDiffs[mateID]));
302 newMateName = mate->GetName() + G4UIcommand::ConvertToString(densityBin);
303 }
304
305 //--- Look if it is the first voxel with this material/densityBin
306 std::pair<G4Material*, G4int> matdens(mate, densityBin);
307
308 auto mppite = newMateDens.find(matdens);
309 if (mppite != newMateDens.cend()) {
310 matInfo* mi = (*mppite).second;
311 mi->fSumdens += dens1;
312 mi->fNvoxels++;
313 fMateIDs[copyNo] = fOriginalMaterials.size() - 1 + mi->fId;
314 // G4cout << copyNo << " mat new again "
315 //<< fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
316 }
317 else {
318 matInfo* mi = new matInfo;
319 mi->fSumdens = dens1;
320 mi->fNvoxels = 1;
321 mi->fId = G4int(newMateDens.size() + 1);
322 newMateDens[matdens] = mi;
323 fMateIDs[copyNo] = fOriginalMaterials.size() - 1 + mi->fId;
324 // G4cout << copyNo << " mat new first "
325 // << fOriginalMaterials.size()-1 + mi->fId << G4endl;
326 }
327 copyNo++;
328 // G4cout << ix << " " << iy << " " << iz
329 // << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
330 // fMateIDs[copyNo] = atoi(cid)-1;
331 }
332 }
333 }
334 }
335
336 //----- Build the list of phantom materials that go to Parameterisation
337 //--- Add original materials
338 for (auto mimite = fOriginalMaterials.cbegin(); mimite != fOriginalMaterials.cend(); ++mimite) {
339 fPhantomMaterials.push_back((*mimite));
340 }
341 //
342 //---- Build and add new materials
343 for (auto mppite = newMateDens.cbegin(); mppite != newMateDens.cend(); ++mppite) {
344 G4double averdens = (*mppite).second->fSumdens / (*mppite).second->fNvoxels;
345 G4double saverdens = G4int(1000.001 * averdens) / 1000.;
346 G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
347 << saverdens << " -> " << G4int(1000 * averdens) << " " << G4int(1000 * averdens) / 1000
348 << " " << G4int(1000 * averdens) / 1000. << G4endl;
349
350 G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
351 << saverdens << " -> " << G4UIcommand::ConvertToString(saverdens) << " Nvoxels "
352 << (*mppite).second->fNvoxels << G4endl;
353 G4String mateName =
354 ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
355 fPhantomMaterials.push_back(
356 BuildMaterialWithChangingDensity((*mppite).first.first, G4float(averdens), mateName));
357 }
358 }
359
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
362 {
363 // build a clinder that encompass the partial phantom built in the example
364 // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
365 //---- Extract number of voxels and voxel dimensions
366
367 G4String contname = "PhantomContainer";
368
369 //----- Define the volume that contains all the voxels
370 G4Tubs* container_tube = new G4Tubs(contname, 0., 200., 100., 0., 360 * deg);
371
372 fContainer_logic = new G4LogicalVolume(container_tube, fAir, contname, 0, 0, 0);
373
374 G4ThreeVector posCentreVoxels(0., 0., 0.);
375 // G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
376 G4RotationMatrix* rotm = new G4RotationMatrix;
377
378 fContainer_phys = new G4PVPlacement(rotm, // rotation
379 posCentreVoxels,
380 fContainer_logic, // The logic volume
381 contname, // Name
382 fWorld_logic, // Mother
383 false, // No op. bool.
384 1); // Copy number
385 }
386
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
388 void DicomPartialDetectorConstruction::ConstructPhantom()
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, fDimX / 2., fDimY / 2., fDimZ / 2.);
399 G4LogicalVolume* phantom_logic = new G4LogicalVolume(phantom_solid, fAir, voxelName, 0, 0, 0);
400 G4bool pVis = 0;
401 if (!pVis) {
402 G4VisAttributes* visAtt = new G4VisAttributes(FALSE);
403 phantom_logic->SetVisAttributes(visAtt);
404 }
405
406 G4double theSmartless = 0.5;
407 // fContainer_logic->SetSmartless( theSmartless );
408 phantom_logic->SetSmartless(theSmartless);
409
410 fPartialPhantomParam = new G4PartialPhantomParameterisation();
411
412 fPartialPhantomParam->SetMaterials(fPhantomMaterials);
413 fPartialPhantomParam->SetVoxelDimensions(fDimX / 2., fDimY / 2., fDimZ / 2.);
414 fPartialPhantomParam->SetNoVoxels(fNoVoxelsX, fNoVoxelsY, fNoVoxelsZ);
415 fPartialPhantomParam->SetMaterialIndices(fMateIDs);
416
417 fPartialPhantomParam->SetFilledIDs(fFilledIDs);
418
419 fPartialPhantomParam->SetFilledMins(fFilledMins);
420
421 fPartialPhantomParam->BuildContainerWalls();
422
423 // G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
424 // G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
425
426 G4PVParameterised* phantom_phys = 0;
427 if (OptimAxis == "kUndefined") {
428 phantom_phys = new G4PVParameterised(voxelName, phantom_logic, fContainer_logic, kUndefined,
429 fNoVoxels, fPartialPhantomParam);
430 }
431 else if (OptimAxis == "kXAxis") {
432 // G4cout << " optim kX " << G4endl;
433 phantom_phys = new G4PVParameterised(voxelName, phantom_logic, fContainer_logic, kXAxis,
434 fNoVoxels, fPartialPhantomParam);
435 fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
436 }
437 else {
438 G4Exception("GmReadPhantomGeometry::ConstructPhantom", "", FatalErrorInArgument,
439 G4String("Wrong argument to parameter \
440 GmReadPhantomGeometry:Phantom:OptimAxis: \
441 Only allowed 'kUndefined' or 'kXAxis', it is: "
442 + OptimAxis)
443 .c_str());
444 }
445
446 phantom_phys->SetRegularStructureId(regStructureID);
447 }
448