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