Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomPhantomZSliceHeader.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 //
26 // $Id: DicomPhantomZSliceHeader.cc 101109 2016-11-07 08:14:53Z gcosmo $
27 //
30 //
31 
32 #include "globals.hh"
33 #include "G4LogicalVolume.hh"
34 #include "G4MaterialTable.hh"
35 #include "G4Material.hh"
36 #include "G4GeometryTolerance.hh"
37 #include "G4NistManager.hh"
38 
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 : fNoVoxelX(0),fNoVoxelY(0),fNoVoxelZ(0),
44  fMinX(0),fMinY(0),fMinZ(0),
45  fMaxX(0),fMaxY(0),fMaxZ(0),
46  fFilename(fname),fSliceLocation(0)
47 {
48 
49 }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 {
54 
55 }
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59  const DicomPhantomZSliceHeader& rhs )
60 {
61  fNoVoxelX = rhs.GetNoVoxelX();
62  fNoVoxelY = rhs.GetNoVoxelY();
63  fNoVoxelZ = rhs.GetNoVoxelZ();
64  fMinX = rhs.GetMinX();
65  fMaxX = rhs.GetMaxX();
66  fMinY = rhs.GetMinY();
67  fMaxY = rhs.GetMaxY();
68  fMinZ = rhs.GetMinZ();
69  fMaxZ = rhs.GetMaxZ();
70  fMaterialNames = rhs.GetMaterialNames();
71  fFilename = rhs.fFilename;
72  fValues = rhs.fValues;
73  fMateIDs = rhs.fMateIDs;
74  fSliceLocation = rhs.fSliceLocation;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 {
80  //----- Read material indices and names
81  G4int nmate;
82  G4String mateindex;
83  G4String matename;
84  fin >> nmate;
85 #ifdef G4VERBOSE
86  G4cout << " DicomPhantomZSliceHeader reading number of materials "
87  << nmate << G4endl;
88 #endif
89 
90  for( G4int im = 0; im < nmate; im++ ){
91  fin >> mateindex >> matename;
92 #ifdef G4VERBOSE
93  //G4cout << " DicomPhantomZSliceHeader reading material " << im << " : "
94  // << mateindex << " " << matename << G4endl;
95 #endif
96 
97  if( ! CheckMaterialExists( matename ) ) {
98  G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader",
99  "A material is found in file that is not built in the C++ code",
100  FatalErrorInArgument, matename.c_str());
101  }
102 
103  fMaterialNames.push_back(matename);
104  }
105 
106  //----- Read number of voxels
107  fin >> fNoVoxelX >> fNoVoxelY >> fNoVoxelZ;
108 #ifdef G4VERBOSE
109  G4cout << " Number of voxels " << fNoVoxelX << " " << fNoVoxelY
110  << " " << fNoVoxelZ << G4endl;
111 #endif
112 
113  //----- Read minimal and maximal extensions (= walls of phantom)
114  fin >> fMinX >> fMaxX;
115  fin >> fMinY >> fMaxY;
116  fin >> fMinZ >> fMaxZ;
117 #ifdef G4VERBOSE
118  G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl
119  << " Extension in Y " << fMinY << " " << fMaxY << G4endl
120  << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
121 #endif
122 
123  fSliceLocation = 0.5*(fMinZ + fMaxZ);
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 G4bool DicomPhantomZSliceHeader::CheckMaterialExists( const G4String& mateName )
128 {
130  std::vector<G4Material*>::const_iterator matite;
131  for( matite = matTab->begin(); matite != matTab->end(); ++matite ) {
132  if( (*matite)->GetName() == mateName ) { return true; }
133  }
134 
136  if( g4mate ) {
137  return false;
138  } else {
139  return true;
140  }
141 
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 {
147  *this = *this + rhs;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152  const DicomPhantomZSliceHeader& rhs )
153 {
154  //----- Check that both slices has the same dimensions
155  if( fNoVoxelX != rhs.GetNoVoxelX()
156  || fNoVoxelY != rhs.GetNoVoxelY() ) {
157  G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
158  !!! Different number of voxels: "
159  << " X= " << fNoVoxelX << " =? " << rhs.GetNoVoxelX()
160  << " Y= " << fNoVoxelY << " =? " << rhs.GetNoVoxelY()
161  << " Z= " << fNoVoxelZ << " =? " << rhs.GetNoVoxelZ()
162  << G4endl;
163  G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader",
164  "",FatalErrorInArgument,"");
165  }
166  //----- Check that both slices has the same extensions
167  if( fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX()
168  || fMinY != rhs.GetMinY() || fMaxY != rhs.GetMaxY() ) {
169  G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
170  !!! Different extensions: "
171  << " Xmin= " << fMinX << " =? " << rhs.GetMinX()
172  << " Xmax= " << fMaxX << " =? " << rhs.GetMaxX()
173  << " Ymin= " << fMinY << " =? " << rhs.GetMinY()
174  << " Ymax= " << fMaxY << " =? " << rhs.GetMaxY()
175  << G4endl;
176  G4Exception("DicomPhantomZSliceHeader::operator+","",
178  }
179 
180  //----- Check that both slices have the same materials
181  std::vector<G4String> fMaterialNames2 = rhs.GetMaterialNames();
182  if( fMaterialNames.size() != fMaterialNames2.size() ) {
183  G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
184  !!! Different number of materials: " << fMaterialNames.size() << " =? "
185  << fMaterialNames2.size() << G4endl;
186  G4Exception("DicomPhantomZSliceHeader::operator+","",
188  }
189  for( unsigned int ii = 0; ii < fMaterialNames.size(); ii++ ) {
190  if( fMaterialNames[ii] != fMaterialNames2[ii] ) {
191  G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
192  !!! Different material number " << ii << " : "
193  << fMaterialNames[ii] << " =? "
194  << fMaterialNames2[ii] << G4endl;
195  G4Exception("DicomPhantomZSliceHeader::operator+","",
197  }
198  }
199 
200  //----- Check that the slices are contiguous in Z
201  if( std::fabs( fMinZ - rhs.GetMaxZ() ) >
203  std::fabs( fMaxZ - rhs.GetMinZ() ) >
205  G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!!\
206  Slices are not contiguous in Z "
207  << " Zmin= " << fMinZ << " & " << rhs.GetMinZ()
208  << " Zmax= " << fMaxZ << " & " << rhs.GetMaxZ()
209  << G4endl;
210  G4Exception("DicomPhantomZSliceHeader::operator+","",
212  }
213 
214  //----- Build slice header copying first one
215  DicomPhantomZSliceHeader temp( *this );
216 
217  //----- Add data from second slice header
218  temp.SetMinZ( std::min( fMinZ, rhs.GetMinZ() ) );
219  temp.SetMaxZ( std::max( fMaxZ, rhs.GetMaxZ() ) );
220  temp.SetNoVoxelZ( fNoVoxelZ + rhs.GetNoVoxelZ() );
221 
222  return temp;
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227 {
228 
229  G4cout << "DicomPhantomZSliceHeader::Dumping Z Slice data to "
230  << fFilename << "..." << G4endl;
231  //sleep(5);
232 
233  // May seen counter-intuitive (dumping to file you are reading from), but
234  // the reason for this is modification slice spacing
235  if(fMateIDs.size() == 0 || fValues.size() == 0) { ReadDataFromFile(); }
236 
237 
238  std::ofstream out;
239  out.open(fFilename.c_str());
240 
241  if(!out) {
242  G4String descript = "DicomPhantomZSliceHeader::DumpToFile: could not open "
243  +fFilename;
244  G4Exception(descript.c_str(),"", FatalException, "");
245  }
246 
247  out << fMaterialNames.size() << std::endl;
248  for(unsigned int i = 0; i < fMaterialNames.size(); ++i) {
249  out << i << " " << fMaterialNames.at(i) << std::endl;
250  }
251 
252  out << fNoVoxelX << " " << fNoVoxelY << " " << fNoVoxelZ << std::endl;
253  out << fMinX << " " << fMaxX << std::endl;
254  out << fMinY << " " << fMaxY << std::endl;
255  out << fMinZ << " " << fMaxZ << std::endl;
256 
257  for(unsigned int i = 0; i < fMateIDs.size(); ++i) { Print(out,fMateIDs.at(i)," "); }
258  for(unsigned int i = 0; i < fValues.size(); ++i) { Print(out,fValues.at(i)," ",6); }
259 
260  out.close();
261 
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 {
267  std::ifstream in;
268  in.open(fFilename.c_str());
269 
270  if(!in) {
271  G4String descript = "DicomPhantomZSliceHeader::DumpToFile: could not open "
272  +fFilename;
273  G4Exception(descript.c_str(),"", FatalException, "");
274  }
275 
276  G4int nMaterials;
277  in >> nMaterials;
278 
279  fMaterialNames.resize(nMaterials,"");
280  for(G4int i = 0; i < nMaterials; ++i) {
281  G4String str1, str2;
282  in >> str1 >> str2;
283  if(!IsInteger(str1)) {
284  G4String descript = "String : " + str1 + " supposed to be integer";
285  G4Exception("DicomPhantomZSliceHeader::ReadDataFromFile - error in \
286  formatting: missing material index","", FatalException,descript.c_str());
287  }
288  G4int index = G4s2n<G4int>(str1);
289  if(index > nMaterials || index < 0) {
290  G4String descript = "Index : " + str1;
291  G4Exception("DicomPhantomZSliceHeader::ReadDataFromFile - error:\
292  bad material index","", FatalException,descript.c_str());
293  }
294  fMaterialNames[index] = str2;
295  }
296 
297  in >> fNoVoxelX >> fNoVoxelY >> fNoVoxelZ;
298 
299  G4double tmpMinX, tmpMinY, tmpMinZ;
300  G4double tmpMaxX, tmpMaxY, tmpMaxZ;
301 
302  in >> tmpMinX >> tmpMaxX;
303  in >> tmpMinY >> tmpMaxY;
304  in >> tmpMinZ >> tmpMaxZ;
305 
306  fMinX = (CheckConsistency(tmpMinX,fMinX,"Min X value")) ?
307  fMinX : ((fMinX == 0) ? tmpMinX : fMinX);
308  fMaxX = (CheckConsistency(tmpMaxX,fMaxX,"Max X value")) ?
309  fMaxX : ((fMaxX == 0) ? tmpMaxX : fMaxX);
310 
311  fMinY = (CheckConsistency(tmpMinY,fMinY,"Min Y value")) ?
312  fMinY : ((fMinY == 0) ? tmpMinY : fMinY);
313  fMaxY = (CheckConsistency(tmpMaxY,fMaxY,"Max Y value")) ?
314  fMaxY : ((fMaxY == 0) ? tmpMaxY : fMaxY);
315 
316  fMinZ = (CheckConsistency(tmpMinZ,fMinZ,"Min Z value")) ?
317  fMinZ : ((fMinZ == 0) ? tmpMinZ : fMinZ);
318  fMaxZ = (CheckConsistency(tmpMaxZ,fMaxZ,"Max Z value")) ?
319  fMaxZ : ((fMaxZ == 0) ? tmpMaxZ : fMaxZ);
320 
321  fMateIDs.clear();
322  fValues.clear();
323  fMateIDs.resize(fNoVoxelY*fNoVoxelZ,std::vector<G4int>(fNoVoxelX,0));
324  fValues.resize(fNoVoxelY*fNoVoxelZ,std::vector<G4double>(fNoVoxelX,0.));
325  for(G4int k = 0; k < fNoVoxelZ; ++k) {
326  for(G4int j = 0; j < fNoVoxelY; ++j) {
327  for(G4int i = 0; i < fNoVoxelX; ++i) {
328  G4int tmpMateID;
329  in >> tmpMateID;
330  G4int row = j*(k+1);
331  fMateIDs[row][i] = tmpMateID;
332  }
333  }
334  }
335 
336  for(G4int k = 0; k < fNoVoxelZ; ++k) {
337  for(G4int j = 0; j < fNoVoxelY; ++j) {
338  for(G4int i = 0; i < fNoVoxelX; ++i) {
339  G4double tmpValue;
340  in >> tmpValue;
341  G4int row = j*(k+1);
342  fValues[row][i] = tmpValue;
343  }
344  }
345  }
346 
347  in.close();
348 }
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
void SetMinZ(const G4double &val)
const std::vector< G4String > & GetMaterialNames() const
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
DicomPhantomZSliceHeader operator+(const DicomPhantomZSliceHeader &rhs)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
DicomPhantomZSliceHeader(const G4String &)
void operator+=(const DicomPhantomZSliceHeader &rhs)
G4GLOB_DLL std::ostream G4cout
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
string fname
Definition: test.py:308
Definition of the DicomPhantomZSliceHeader class.
double G4double
Definition: G4Types.hh:76
void SetMaxZ(const G4double &val)
void SetNoVoxelZ(const G4int &val)
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr