Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomPartialDetectorConstruction.cc
Go to the documentation of this file.
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 //
28 //
29 
30 #include "globals.hh"
31 
32 #include "DicomPartialDetectorConstruction.hh"
33 
34 #include "G4Box.hh"
35 #include "G4LogicalVolume.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4PVPlacement.hh"
38 #include "G4PVParameterised.hh"
39 #include "G4Material.hh"
40 #include "G4Element.hh"
41 #include "G4VisAttributes.hh"
42 #include "G4Colour.hh"
43 #include "G4ios.hh"
44 #include "G4PartialPhantomParameterisation.hh"
45 #include "G4NistManager.hh"
46 #include "G4UIcommand.hh"
47 #include "G4Tubs.hh"
48 
49 #include "G4tgbVolumeMgr.hh"
50 #include "G4tgbMaterialMgr.hh"
51 #include "G4tgrMessenger.hh"
52 #include "G4tgrMessenger.hh"
53 
54 #include "G4SystemOfUnits.hh"
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 DicomPartialDetectorConstruction::DicomPartialDetectorConstruction()
58 {
59 
60  fPartialPhantomParam = 0;
61 
62  fNVoxels = 0;
63  fDimX = 0;
64  fDimY = 0;
65  fDimZ = 0;
66  fOffsetX = 0;
67  fOffsetY = 0;
68  fOffsetZ = 0;
69  fMateIDs = 0;
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction()
74 {
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 G4VPhysicalVolume* DicomPartialDetectorConstruction::Construct()
79 {
80  InitialisationOfMaterials();
81 
82  //----- Build world
83  G4double worldXDimension = 1.*m;
84  G4double worldYDimension = 1.*m;
85  G4double worldZDimension = 1.*m;
86 
87  G4Box* world_box = new G4Box( "WorldSolid",
88  worldXDimension,
89  worldYDimension,
90  worldZDimension );
91 
92  fWorld_logic = new G4LogicalVolume( world_box,
93  fAir,
94  "WorldLogical",
95  0, 0, 0 );
96 
97  G4VPhysicalVolume* world_phys = new G4PVPlacement( 0,
98  G4ThreeVector(0,0,0),
99  "World",
100  fWorld_logic,
101  0,
102  false,
103  0 );
104 
105  ReadPhantomData();
106 
107  ConstructPhantomContainer();
108  ConstructPhantom();
109 
110  return world_phys;
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 void DicomPartialDetectorConstruction::ReadPhantomData()
115 {
116  std::ifstream finDF("Data.dat");
117  G4String fname;
118  if(finDF.good() != 1 ) {
119  G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData",
120  "",
121  FatalErrorInArgument,
122  " Problem reading data file: Data.dat");
123  }
124 
125  G4int compression;
126  finDF >> compression; // not used here
127 
128  G4int numFiles;
129  finDF >> numFiles; // only 1 file supported for the moment
130  if( numFiles != 1 ) {
131  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
132  "",
133  FatalErrorInArgument,
134  "More than 1 DICOM file is not supported");
135  }
136  for(G4int i = 0; i < numFiles; i++ ) {
137  finDF >> fname;
138  //--- Read one data file
139  fname += ".g4pdcm";
140  ReadPhantomDataFile(fname);
141  }
142 
143  finDF.close();
144 
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
148 void DicomPartialDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
149 {
150  // G4String filename = "phantom.g4pdcm";
151 #ifdef G4VERBOSE
152  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
153 #endif
154  std::ifstream fin(fname.c_str(), std::ios_base::in);
155  if( !fin.is_open() ) {
156  G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile",
157  "",
158  FatalException,
159  G4String("File not found " + fname).c_str());
160  }
161  G4int nMaterials;
162  fin >> nMaterials;
163  G4String stemp;
164  G4int nmate;
165  for( G4int ii = 0; ii < nMaterials; ii++ ){
166  fin >> nmate >> stemp;
167  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate " << ii << " = " << nmate << " mate " << stemp << G4endl;
168  if( ii != nmate ) G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
169  "",
170  FatalErrorInArgument,
171  "Material number should be in increasing order: wrong material number ");
172 
173  }
174 
175  fin >> fNVoxelX >> fNVoxelY >> fNVoxelZ;
176  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z " << fNVoxelX << " " << fNVoxelY << " " << fNVoxelZ << G4endl;
177  fin >> fOffsetX >> fDimX;
178  fDimX = (fDimX - fOffsetX)/fNVoxelX;
179  fin >> fOffsetY >> fDimY;
180  fDimY = (fDimY - fOffsetY)/fNVoxelY;
181  fin >> fOffsetZ >> fDimZ;
182  fDimZ = (fDimZ - fOffsetZ)/fNVoxelZ;
183  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimX " << fDimX << " fOffsetX " << fOffsetX << G4endl;
184  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY " << fDimY << " fOffsetY " << fOffsetY << G4endl;
185  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ " << fDimZ << " fOffsetZ " << fOffsetZ << G4endl;
186 
187  //--- Read voxels that are filled
188  fNVoxels = 0;
189  // G4bool* isFilled = new G4bool[fNVoxelX*fNVoxelY*fNVoxelZ];
190  // fFilledIDs = new size_t[fNVoxelZ*fNVoxelY+1];
191  //? fFilledIDs.insert(0);
192  int ifxmin1, ifxmax1;
193  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
194  std::map< G4int, G4int > ifmin, ifmax;
195  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
196  fin >> ifxmin1 >> ifxmax1;
197  // check coherence ...
198 
199  ifmin[iy] = ifxmin1;
200  ifmax[iy] = ifxmax1;
201  G4int ifxdiff = ifxmax1-ifxmin1+1;
202  if( ifxmax1 == -1 && ifxmin1 == -1 ) ifxdiff = 0;
203  // fFilledIDs.insert(std::multimap<G4int,size_t>::value_type(ifxdiff+fNVoxels, ifxmin1));
204  fFilledIDs.insert(std::pair<size_t,size_t>(ifxdiff+fNVoxels-1, ifxmin1));
205  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert FilledIDs " << ifxdiff+fNVoxels-1 << " min " << ifxmin1 << " N= " << fFilledIDs.size() << G4endl;
206  //filledIDs[iz*fNVoxelY+iy+1] = ifxmax1-ifxmin1+1;
207  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
208  // G4cout << "REGNAV DicomPartialDetectorConstruction::ReadPhantomData checking to add voxel iz " << iz << " iy " << iy << " ix " << ix << " min= " << ifxmin1 << " max= " << ifxmax1 << G4endl;
209  // G4int copyNo = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
210  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
211  // isFilled[copyNo] = true;
212  fNVoxels++;
213  } else {
214  // isFilled[copyNo] = false;
215  }
216  }
217  }
218  fFilledMins[iz] = ifmin;
219  fFilledMaxs[iz] = ifmax;
220  }
221 
222  //--- Read material IDs
223  G4int mateID1;
224  fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ];
225  G4int copyNo = 0;
226  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
227  std::map< G4int, G4int > ifmin = fFilledMins[iz];
228  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
229  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
230  ifxmin1 = ifmin[iy];
231  ifxmax1 = ifmax[iy];
232  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
233  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
234  fin >> mateID1;
235  if( mateID1 < 0 || mateID1 >= nMaterials ) {
236  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
237  "",
238  FatalException,
239  G4String("Wrong index in phantom file: It should be between 0 and "
240  + G4UIcommand::ConvertToString(nMaterials-1)
241  + ", while it is "
242  + G4UIcommand::ConvertToString(mateID1)).c_str());
243  }
244  fMateIDs[copyNo] = mateID1;
245  copyNo++;
246  }
247  }
248  }
249  }
250 
251  ReadVoxelDensitiesPartial( fin );
252 
253  fin.close();
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
257 void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial( std::ifstream& fin )
258 {
259  std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
260  std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
261  for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
262  densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
263  }
264 
265  //----- Define density differences (maximum density difference to create a new material)
266  char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
267  G4double densityDiff = -1.;
268  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
269 
270  std::map<G4int,G4double> densityDiffs;
271  if( densityDiff != -1. ) {
272  for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
273  densityDiffs[ii] = densityDiff; //currently all materials with same step
274  }
275  }else {
276  if( fMaterials.size() == 0 ) { // do it only for first slice
277  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
278  fMaterials.push_back( fOriginalMaterials[ii] );
279  }
280  }
281  }
282 
283  // densityDiffs[0] = 0.0001; //fAir
284 
285  //--- Calculate the average material density for each material/density bin
286  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
287 
288  G4double dens1;
289  size_t ifxmin1, ifxmax1;
290 
291  //---- Read the material densities
292  G4int copyNo = 0;
293  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
294  std::map< G4int, G4int > ifmin = fFilledMins[iz];
295  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
296  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
297  ifxmin1 = ifmin[iy];
298  ifxmax1 = ifmax[iy];
299  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
300  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
301  fin >> dens1;
302  // G4cout << ix << " " << iy << " " << iz << " filling fMateIDs " << copyNo << " = " << atoi(stemp.c_str())-1 << " " << stemp << G4endl;
303 
304  //--- store the minimum and maximum density for each material (just for printing)
305  mpite = densiMinMax.find( fMateIDs[copyNo] );
306  if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
307  if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
308 
309  //--- Get material from original list of material in file
310  int mateID = fMateIDs[copyNo];
311  G4Material* mate = fOriginalMaterials[mateID];
312  // G4cout << copyNo << " mateID " << mateID << G4endl;
313  //--- Check if density is equal to the original material density
314  if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 ) continue;
315 
316  //--- Build material name with fOriginalMaterials name + density
317  // float densityBin = densityDiffs[mateID] * (G4int(dens1/densityDiffs[mateID])+0.5);
318  G4String newMateName = mate->GetName();
319 
320  G4int densityBin = 0;
321  // G4cout << " densityBin " << densityBin << " " << dens1 << " " <<densityDiffs[mateID] << G4endl;
322 
323  if( densityDiff != -1.) {
324  densityBin = (G4int(dens1/densityDiffs[mateID]));
325  newMateName = mate->GetName()+G4UIcommand::ConvertToString(densityBin);
326  }
327 
328  //--- Look if it is the first voxel with this material/densityBin
329  std::pair<G4Material*,G4int> matdens(mate, densityBin );
330 
331  std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite = newMateDens.find( matdens );
332  if( mppite != newMateDens.end() ){
333  matInfo* mi = (*mppite).second;
334  mi->sumdens += dens1;
335  mi->nvoxels++;
336  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
337  // G4cout << copyNo << " mat new again " << fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
338  } else {
339  matInfo* mi = new matInfo;
340  mi->sumdens = dens1;
341  mi->nvoxels = 1;
342  mi->id = newMateDens.size()+1;
343  newMateDens[matdens] = mi;
344  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
345  // G4cout << copyNo << " mat new first " << fOriginalMaterials.size()-1 + mi->id << G4endl;
346  }
347  copyNo++;
348  // G4cout << ix << " " << iy << " " << iz << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
349  // fMateIDs[copyNo] = atoi(cid)-1;
350  }
351  }
352  }
353  }
354 
355  //----- Build the list of phantom materials that go to Parameterisation
356  //--- Add original materials
357  std::vector<G4Material*>::const_iterator mimite;
358  for( mimite = fOriginalMaterials.begin(); mimite != fOriginalMaterials.end(); mimite++ ){
359  fPhantomMaterials.push_back( (*mimite) );
360  }
361  //
362  //---- Build and add new materials
363  std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite;
364  for( mppite= newMateDens.begin(); mppite != newMateDens.end(); mppite++ ){
365  G4double averdens = (*mppite).second->sumdens/(*mppite).second->nvoxels;
366  G4double saverdens = G4int(1000.001*averdens)/1000.;
367  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> " << saverdens << " -> " << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000 << " " << G4int(1000*averdens)/1000. << G4endl;
368 
369  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> " << saverdens << " -> " << G4UIcommand::ConvertToString(saverdens) << " Nvoxels " <<(*mppite).second->nvoxels << G4endl;
370  G4String mateName = ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
371  fPhantomMaterials.push_back( BuildMaterialWithChangingDensity( (*mppite).first.first, averdens, mateName ) );
372  }
373 
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
378 {
379  //build a clinder that encompass the partial phantom built in the example
380  // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
381  //---- Extract number of voxels and voxel dimensions
382 
383  G4String contname= "PhantomContainer";
384 
385  //----- Define the volume that contains all the voxels
386  G4Tubs* container_tube = new G4Tubs(contname,0., 200., 100., 0., 360*deg);
387 
388  fContainer_logic =
389  new G4LogicalVolume( container_tube,
390  fAir,
391  contname,
392  0, 0, 0 );
393 
394  G4ThreeVector posCentreVoxels(0.,0.,0.);
395  //G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
396  G4RotationMatrix* rotm = new G4RotationMatrix;
397 
398  fContainer_phys =
399  new G4PVPlacement(rotm, // rotation
400  posCentreVoxels,
401  fContainer_logic, // The logic volume
402  contname, // Name
403  fWorld_logic, // Mother
404  false, // No op. bool.
405  1); // Copy number
406 
407 }
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410 void DicomPartialDetectorConstruction::ConstructPhantom()
411 {
412  G4String OptimAxis = "kXAxis";
413  G4bool bSkipEqualMaterials = 0;
414  G4int regStructureID = 2;
415 
416  G4ThreeVector posCentreVoxels(0.,0.,0.);
417 
418  //----- Build phantom
419  G4String voxelName = "phantom";
420  G4Box* phantom_solid = new G4Box(voxelName, fDimX/2., fDimY/2., fDimZ/2. );
421  G4LogicalVolume* phantom_logic =
422  new G4LogicalVolume( phantom_solid,
423  fAir,
424  voxelName,
425  0, 0, 0 );
426  G4bool pVis = 0;
427  if( !pVis ) {
428  G4VisAttributes* visAtt = new G4VisAttributes( FALSE );
429  phantom_logic->SetVisAttributes( visAtt );
430  }
431 
432  G4double theSmartless = 0.5;
433  // fContainer_logic->SetSmartless( theSmartless );
434  phantom_logic->SetSmartless( theSmartless );
435 
436  fPartialPhantomParam = new G4PartialPhantomParameterisation();
437 
438  fPartialPhantomParam->SetMaterials( fPhantomMaterials );
439  fPartialPhantomParam->SetVoxelDimensions( fDimX/2., fDimY/2., fDimZ/2. );
440  fPartialPhantomParam->SetNoVoxel( fNVoxelX, fNVoxelY, fNVoxelZ );
441  fPartialPhantomParam->SetMaterialIndices( fMateIDs );
442 
443  fPartialPhantomParam->SetFilledIDs(fFilledIDs);
444 
445  fPartialPhantomParam->SetFilledMins(fFilledMins);
446 
447  fPartialPhantomParam->BuildContainerWalls();
448 
449  // G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
450  // G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
451 
452  G4PVParameterised * phantom_phys = 0;
453  if( OptimAxis == "kUndefined" ) {
454  phantom_phys = new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
455  kUndefined, fNVoxels, fPartialPhantomParam);
456  } else if( OptimAxis == "kXAxis" ) {
457  // G4cout << " optim kX " << G4endl;
458  phantom_phys = new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
459  // kXAxis, fNVoxelX*fNVoxelY*fNVoxelZ, fPartialPhantomParam);
460  kXAxis, fNVoxels, fPartialPhantomParam);
461  fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
462  } else {
463  G4Exception("GmReadPhantomGeometry::ConstructPhantom",
464  "",
465  FatalErrorInArgument,
466  (G4String("Wrong argument to parameter GmReadPhantomGeometry:Phantom:OptimAxis: Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str()));
467  }
468 
469  phantom_phys->SetRegularStructureId(regStructureID);
470 
471  // G4cout << " GmReadPhantomGeometry phantom_phys " << phantom_phys << " trans " << phantom_phys->GetTranslation() << G4endl;
472  //- fPartialPhantomParam->BuildContainerSolid(cont_phys);
473 
474  // G4cout << " GmReadPhantomGeometry::constructPhantom ended " << G4endl;
475 
476 }
Generated on Sat May 25 2013 14:32:57 for Geant4 by   doxygen 1.8.4