| 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/DicomIntersectVolu
28 /// \brief Implementation of the DicomIntersec
29 //
30
31 #include "DicomIntersectVolume.hh"
32
33 #include "G4LogicalVolume.hh"
34 #include "G4LogicalVolumeStore.hh"
35 #include "G4Material.hh"
36 #include "G4PVParameterised.hh"
37 #include "G4PhantomParameterisation.hh"
38 #include "G4PhysicalVolumeStore.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UIcmdWithAString.hh"
41 #include "G4UIcommand.hh"
42 #include "G4VPhysicalVolume.hh"
43 #include "G4VSolid.hh"
44 #include "G4tgbVolume.hh"
45 #include "G4tgrSolid.hh"
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48 DicomIntersectVolume::DicomIntersectVolume()
49 : G4UImessenger(), fG4VolumeCmd(0), fSolid(0
50 {
51 fUserVolumeCmd = new G4UIcmdWithAString("/di
52 fUserVolumeCmd->SetGuidance(
53 "Intersects a phantom with a user-defined
54 "and outputs the voxels that are totally i
55 " a new phantom file. It must have the par
56 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_
57 fUserVolumeCmd->SetParameterName("choice", t
58 fUserVolumeCmd->AvailableForStates(G4State_I
59
60 fG4VolumeCmd = new G4UIcmdWithAString("/dico
61 fG4VolumeCmd->SetGuidance(
62 "Intersects a phantom with a user-defined
63 " and outputs the voxels that are totally
64 " a new phantom file. It must have the par
65 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_
66 fG4VolumeCmd->SetParameterName("choice", tru
67 fG4VolumeCmd->AvailableForStates(G4State_Idl
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71 DicomIntersectVolume::~DicomIntersectVolume()
72 {
73 if (fUserVolumeCmd) delete fUserVolumeCmd;
74 if (fG4VolumeCmd) delete fG4VolumeCmd;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oo
78 void DicomIntersectVolume::SetNewValue(G4UIcom
79 {
80 G4AffineTransform theVolumeTransform;
81
82 if (command == fUserVolumeCmd) {
83 std::vector<G4String> params = GetWordsInS
84 if (params.size() < 8) {
85 G4Exception("DicomIntersectVolume::SetNe
86 " There must be at least 8 p
87 "ANG_X ANG_Y ANG_Z SOLID_PAR
88 FatalErrorInArgument,
89 G4String("Number of paramete
90 + G4UIcommand::Conv
91 .c_str());
92 }
93
94 //----- Build G4VSolid
95 BuildUserSolid(params);
96
97 //----- Calculate volume inverse 3D transf
98 G4ThreeVector pos = G4ThreeVector(G4UIcomm
99 G4UIcomm
100 G4UIcomm
101 G4RotationMatrix* rotmat = new G4RotationM
102 std::vector<G4double> angles;
103 rotmat->rotateX(G4UIcommand::ConvertToDoub
104 rotmat->rotateY(G4UIcommand::ConvertToDoub
105 rotmat->rotateY(G4UIcommand::ConvertToDoub
106 theVolumeTransform = G4AffineTransform(rot
107 }
108 else if (command == fG4VolumeCmd) {
109 std::vector<G4String> params = GetWordsInS
110 if (params.size() != 1)
111 G4Exception("DicomIntersectVolume::SetNe
112 G4String("Command: " + comma
113 + " " + newValues +
114 .c_str());
115
116 //----- Build G4VSolid
117 BuildG4Solid(params);
118
119 //----- Calculate volume inverse 3D transf
120 G4VPhysicalVolume* pv = GetPhysicalVolumes
121
122 theVolumeTransform = G4AffineTransform(pv-
123 }
124
125 //----- Calculate relative phantom - volume
126 G4PhantomParameterisation* thePhantomParam =
127
128 G4RotationMatrix* rotph = new G4RotationMatr
129 // assumes the phantom mother is not rotated
130 G4AffineTransform thePhantomTransform(rotph,
131 // assumes the phantom mother is not transla
132
133 G4AffineTransform theTransform = theVolumeTr
134
135 G4ThreeVector axisX(1., 0., 0.);
136 G4ThreeVector axisY(0., 1., 0.);
137 G4ThreeVector axisZ(0., 0., 1.);
138 theTransform.ApplyAxisTransform(axisX);
139 theTransform.ApplyAxisTransform(axisY);
140 theTransform.ApplyAxisTransform(axisZ);
141
142 //----- Write phantom header
143 G4String thePhantomFileName = "phantom.g4pdc
144 fout.open(thePhantomFileName);
145 std::vector<G4Material*> materials = thePhan
146 fout << materials.size() << G4endl;
147 for (unsigned int ii = 0; ii < materials.siz
148 fout << ii << " " << materials[ii]->GetNam
149 }
150
151 //----- Loop to pantom voxels
152 G4int nx = G4int(thePhantomParam->GetNoVoxel
153 G4int ny = G4int(thePhantomParam->GetNoVoxel
154 G4int nz = G4int(thePhantomParam->GetNoVoxel
155 G4int nxy = nx * ny;
156 fVoxelIsInside = new G4bool[nx * ny * nz];
157 G4double voxelHalfWidthX = thePhantomParam->
158 G4double voxelHalfWidthY = thePhantomParam->
159 G4double voxelHalfWidthZ = thePhantomParam->
160
161 //----- Write phantom dimensions and limits
162 fout << nx << " " << ny << " " << nz << G4en
163 fout << -voxelHalfWidthX * nx + thePhantomTr
164 << voxelHalfWidthX * nx + thePhantomTra
165 fout << -voxelHalfWidthY * ny + thePhantomTr
166 << voxelHalfWidthY * ny + thePhantomTra
167 fout << -voxelHalfWidthZ * nz + thePhantomTr
168 << voxelHalfWidthZ * nz + thePhantomTra
169
170 for (G4int iz = 0; iz < nz; ++iz) {
171 for (G4int iy = 0; iy < ny; ++iy) {
172 G4bool bPrevVoxelInside = true;
173 G4bool b1VoxelFoundInside = false;
174 G4int firstVoxel = -1;
175 G4int lastVoxel = -1;
176 for (G4int ix = 0; ix < nx; ++ix) {
177 G4ThreeVector voxelCentre((-nx + ix *
178 (-ny + iy *
179 (-nz + iz *
180 theTransform.ApplyPointTransform(voxel
181 G4bool bVoxelIsInside = true;
182 for (G4int ivx = -1; ivx <= 1; ivx +=
183 for (G4int ivy = -1; ivy <= 1; ivy +
184 for (G4int ivz = -1; ivz <= 1; ivz
185 G4ThreeVector voxelPoint = voxel
186 + ivy
187 + ivz
188 if (fSolid->Inside(voxelPoint) =
189 bVoxelIsInside = false;
190 break;
191 }
192 else {
193 }
194 }
195 if (!bVoxelIsInside) break;
196 }
197 if (!bVoxelIsInside) break;
198 }
199
200 G4int copyNo = ix + nx * iy + nxy * iz
201 if (bVoxelIsInside) {
202 if (!bPrevVoxelInside) {
203 G4Exception("DicomIntersectVolume:
204 "Volume cannot interse
205 "please use other voxe
206 }
207 if (!b1VoxelFoundInside) {
208 firstVoxel = ix;
209 b1VoxelFoundInside = true;
210 }
211 lastVoxel = ix;
212 fVoxelIsInside[copyNo] = true;
213 }
214 else {
215 fVoxelIsInside[copyNo] = false;
216 }
217 }
218 fout << firstVoxel << " " << lastVoxel <
219 }
220 }
221
222 //----- Now write the materials
223 for (G4int iz = 0; iz < nz; ++iz) {
224 for (G4int iy = 0; iy < ny; ++iy) {
225 G4bool b1xFound = false;
226 for (G4int ix = 0; ix < nx; ++ix) {
227 size_t copyNo = ix + ny * iy + nxy * i
228 if (fVoxelIsInside[copyNo]) {
229 fout << thePhantomParam->GetMaterial
230 b1xFound = true;
231 }
232 }
233 if (b1xFound) fout << G4endl;
234 }
235 }
236
237 // Now write densities
238 for (G4int iz = 0; iz < nz; ++iz) {
239 for (G4int iy = 0; iy < ny; ++iy) {
240 G4bool b1xFound = false;
241 for (G4int ix = 0; ix < nx; ++ix) {
242 size_t copyNo = ix + ny * iy + nxy * i
243 if (fVoxelIsInside[copyNo]) {
244 fout << thePhantomParam->GetMaterial
245 b1xFound = true;
246 }
247 }
248 if (b1xFound) fout << G4endl;
249 }
250 }
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oo
254 void DicomIntersectVolume::BuildUserSolid(std:
255 {
256 for (G4int ii = 0; ii < 6; ++ii)
257 params.erase(params.begin());
258 // take otu position and rotation angles
259 params.insert(params.begin(), ":SOLID");
260 params.insert(params.begin(), params[1]);
261 G4tgrSolid* tgrSolid = new G4tgrSolid(params
262 G4tgbVolume* tgbVolume = new G4tgbVolume();
263 fSolid = tgbVolume->FindOrConstructG4Solid(t
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oo
267 void DicomIntersectVolume::BuildG4Solid(std::v
268 {
269 fSolid = GetLogicalVolumes(params[0], 1, 1)[
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oo
273 G4PhantomParameterisation* DicomIntersectVolum
274 {
275 G4PhantomParameterisation* paramreg = 0;
276
277 G4PhysicalVolumeStore* pvs = G4PhysicalVolum
278 for (auto cite = pvs->cbegin(); cite != pvs-
279 if (IsPhantomVolume(*cite)) {
280 const G4PVParameterised* pvparam = stati
281 G4VPVParameterisation* param = pvparam->
282 // if( static_cast<const G4PhantomPar
283 // if( static_cast<const G4PhantomPar
284 // G4cout << "G4PhantomParameterisa
285 // << (*cite)->GetName() << G4endl;
286 paramreg = static_cast<G4PhantomParamete
287 }
288 }
289
290 if (!paramreg && bMustExist)
291 G4Exception("DicomIntersectVolume::GetPhan
292 " No G4PhantomParameterisation
293
294 return paramreg;
295 }
296
297 //....oooOO0OOooo........oooOO0OOooo........oo
298 std::vector<G4VPhysicalVolume*> DicomIntersect
299
300 {
301 std::vector<G4VPhysicalVolume*> vvolu;
302 std::string::size_type ial = name.rfind(":")
303 G4String volname = "";
304 G4int volcopy = 0;
305 if (ial != std::string::npos) {
306 std::string::size_type ial2 = name.rfind("
307 if (ial2 != std::string::npos) {
308 G4Exception("DicomIntersectVolume::GetPh
309 G4String("Name corresponds t
310 }
311 else {
312 volname = name.substr(0, ial);
313 volcopy = G4UIcommand::ConvertToInt(name
314 }
315 }
316 else {
317 volname = name;
318 volcopy = -1;
319 }
320
321 G4PhysicalVolumeStore* pvs = G4PhysicalVolum
322 for (auto citepv = pvs->cbegin(); citepv !=
323 if (volname == (*citepv)->GetName() && ((*
324 vvolu.push_back(*citepv);
325 }
326 }
327
328 //----- Check that at least one volume was f
329 if (vvolu.size() == 0) {
330 if (exists) {
331 G4Exception(" DicomIntersectVolume::GetP
332 G4String("No physical volume
333 }
334 else {
335 G4cerr << "!!WARNING: DicomIntersectVolu
336 << " no physical volume found wit
337 }
338 }
339
340 if (nVols != -1 && G4int(vvolu.size()) != nV
341 G4Exception("DicomIntersectVolume::GetLogi
342 "Wrong number of physical volu
343 ("Number of physical volumes "
344 + ", requesting " + G4UIcomma
345 .c_str());
346 }
347
348 return vvolu;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352 G4bool DicomIntersectVolume::IsPhantomVolume(G
353 {
354 EAxis axis;
355 G4int nReplicas;
356 G4double width, offset;
357 G4bool consuming;
358 pv->GetReplicationData(axis, nReplicas, widt
359 EVolume type = (consuming) ? kReplica : kPar
360 if (type == kParameterised && pv->GetRegular
361 return TRUE;
362 }
363 else {
364 return FALSE;
365 }
366 }
367
368 //....oooOO0OOooo........oooOO0OOooo........oo
369 std::vector<G4LogicalVolume*> DicomIntersectVo
370
371 {
372 // G4cout << "GetLogicalVolumes " << name <
373 std::vector<G4LogicalVolume*> vvolu;
374 G4int ial = G4int(name.rfind(":"));
375 if (ial != -1) {
376 G4Exception("DicomIntersectVolume::GetLogi
377 G4String("Name corresponds to
378 }
379
380 G4LogicalVolumeStore* lvs = G4LogicalVolumeS
381 for (auto citelv = lvs->cbegin(); citelv !=
382 if (name == (*citelv)->GetName()) {
383 vvolu.push_back(*citelv);
384 }
385 }
386
387 //----- Check that at least one volume was f
388 if (vvolu.size() == 0) {
389 if (exists) {
390 G4Exception("DicomIntersectVolume::GetLo
391 ("no logical volume found wi
392 }
393 else {
394 G4Exception("DicomIntersectVolume::GetLo
395 ("no logical volume found w
396 }
397 }
398
399 if (nVols != -1 && G4int(vvolu.size()) != nV
400 G4Exception("DicomIntersectVolume::GetLogi
401 FatalErrorInArgument,
402 ("Number of logical volumes "
403 + ", requesting " + G4UIcomma
404 .c_str());
405 }
406
407 return vvolu;
408 }
409
410 //....oooOO0OOooo........oooOO0OOooo........oo
411 std::vector<G4String> DicomIntersectVolume::Ge
412 {
413 std::vector<G4String> wordlist;
414
415 //---------- Read a line of file:
416 //----- Clear wordlist
417 G4int ii;
418 const char* cstr = stemp.c_str();
419 G4int siz = G4int(stemp.length());
420 G4int istart = 0;
421 G4int nQuotes = 0;
422 G4bool lastIsBlank = false;
423 G4bool lastIsQuote = false;
424 for (ii = 0; ii < siz; ++ii) {
425 if (cstr[ii] == '\"') {
426 if (lastIsQuote) {
427 G4Exception("GmGenUtils:GetWordsFromSt
428 ("There cannot be two quot
429 }
430 if (nQuotes % 2 == 1) {
431 // close word
432 wordlist.push_back(stemp.substr(istart
433 // G4cout << "GetWordsInString
434 //<< wordlist[wordlist.size()-1] << "
435 //<< " ii " << ii << G4endl;
436 }
437 else {
438 istart = ii + 1;
439 }
440 ++nQuotes;
441 lastIsQuote = true;
442 lastIsBlank = false;
443 }
444 else if (cstr[ii] == ' ') {
445 // G4cout << "GetWordsInStrin
446 //<< nQuotes << " lastIsBlank " << lastI
447 if (nQuotes % 2 == 0) {
448 if (!lastIsBlank && !lastIsQuote) {
449 wordlist.push_back(stemp.substr(ista
450 // G4cout << "GetWordsInStr
451 //<< wordlist[wordlist.size()-1] <<
452 //<< lastIsBlank << G4endl;
453 }
454
455 istart = ii + 1;
456 lastIsQuote = false;
457 lastIsBlank = true;
458 }
459 }
460 else {
461 if (ii == siz - 1) {
462 wordlist.push_back(stemp.substr(istart
463 // G4cout << "GetWordsI
464 //<< wordlist[wordlist.size()-1] << "
465 }
466 lastIsQuote = false;
467 lastIsBlank = false;
468 }
469 }
470 if (nQuotes % 2 == 1) {
471 G4Exception("GmGenUtils:GetWordsFromString
472 ("unbalanced quotes in line "
473 }
474
475 return wordlist;
476 }
477