| 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