| 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 medical/DICOM/src/DicomIntersectVolume.cc
28 /// \brief Implementation of the DicomIntersectVolume class
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........oooOO0OOooo........oooOO0OOooo......
48 DicomIntersectVolume::DicomIntersectVolume()
49 : G4UImessenger(), fG4VolumeCmd(0), fSolid(0), fPhantomMinusCorner(), fVoxelIsInside(0)
50 {
51 fUserVolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume", this);
52 fUserVolumeCmd->SetGuidance(
53 "Intersects a phantom with a user-defined volume "
54 "and outputs the voxels that are totally inside the intersection as"
55 " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
56 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
57 fUserVolumeCmd->SetParameterName("choice", true);
58 fUserVolumeCmd->AvailableForStates(G4State_Idle);
59
60 fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume", this);
61 fG4VolumeCmd->SetGuidance(
62 "Intersects a phantom with a user-defined volume"
63 " and outputs the voxels that are totally inside the intersection as "
64 " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
65 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
66 fG4VolumeCmd->SetParameterName("choice", true);
67 fG4VolumeCmd->AvailableForStates(G4State_Idle);
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 DicomIntersectVolume::~DicomIntersectVolume()
72 {
73 if (fUserVolumeCmd) delete fUserVolumeCmd;
74 if (fG4VolumeCmd) delete fG4VolumeCmd;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 void DicomIntersectVolume::SetNewValue(G4UIcommand* command, G4String newValues)
79 {
80 G4AffineTransform theVolumeTransform;
81
82 if (command == fUserVolumeCmd) {
83 std::vector<G4String> params = GetWordsInString(newValues);
84 if (params.size() < 8) {
85 G4Exception("DicomIntersectVolume::SetNewValue",
86 " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
87 "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
88 FatalErrorInArgument,
89 G4String("Number of parameters given = "
90 + G4UIcommand::ConvertToString(G4int(params.size())))
91 .c_str());
92 }
93
94 //----- Build G4VSolid
95 BuildUserSolid(params);
96
97 //----- Calculate volume inverse 3D transform
98 G4ThreeVector pos = G4ThreeVector(G4UIcommand::ConvertToDouble(params[0]),
99 G4UIcommand::ConvertToDouble(params[1]),
100 G4UIcommand::ConvertToDouble(params[2]));
101 G4RotationMatrix* rotmat = new G4RotationMatrix;
102 std::vector<G4double> angles;
103 rotmat->rotateX(G4UIcommand::ConvertToDouble(params[3]));
104 rotmat->rotateY(G4UIcommand::ConvertToDouble(params[4]));
105 rotmat->rotateY(G4UIcommand::ConvertToDouble(params[5]));
106 theVolumeTransform = G4AffineTransform(rotmat, pos).Invert();
107 }
108 else if (command == fG4VolumeCmd) {
109 std::vector<G4String> params = GetWordsInString(newValues);
110 if (params.size() != 1)
111 G4Exception("DicomIntersectVolume::SetNewValue", "", FatalErrorInArgument,
112 G4String("Command: " + command->GetCommandPath() + "/" + command->GetCommandName()
113 + " " + newValues + " needs 1 argument: VOLUME_NAME")
114 .c_str());
115
116 //----- Build G4VSolid
117 BuildG4Solid(params);
118
119 //----- Calculate volume inverse 3D transform
120 G4VPhysicalVolume* pv = GetPhysicalVolumes(params[0], 1, 1)[0];
121
122 theVolumeTransform = G4AffineTransform(pv->GetFrameRotation(), pv->GetFrameTranslation());
123 }
124
125 //----- Calculate relative phantom - volume 3D transform
126 G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
127
128 G4RotationMatrix* rotph = new G4RotationMatrix();
129 // assumes the phantom mother is not rotated !!!
130 G4AffineTransform thePhantomTransform(rotph, G4ThreeVector());
131 // assumes the phantom mother is not translated !!!
132
133 G4AffineTransform theTransform = theVolumeTransform * thePhantomTransform;
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.g4pdcm";
144 fout.open(thePhantomFileName);
145 std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
146 fout << materials.size() << G4endl;
147 for (unsigned int ii = 0; ii < materials.size(); ++ii) {
148 fout << ii << " " << materials[ii]->GetName() << G4endl;
149 }
150
151 //----- Loop to pantom voxels
152 G4int nx = G4int(thePhantomParam->GetNoVoxelsX());
153 G4int ny = G4int(thePhantomParam->GetNoVoxelsY());
154 G4int nz = G4int(thePhantomParam->GetNoVoxelsZ());
155 G4int nxy = nx * ny;
156 fVoxelIsInside = new G4bool[nx * ny * nz];
157 G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
158 G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
159 G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
160
161 //----- Write phantom dimensions and limits
162 fout << nx << " " << ny << " " << nz << G4endl;
163 fout << -voxelHalfWidthX * nx + thePhantomTransform.NetTranslation().x() << " "
164 << voxelHalfWidthX * nx + thePhantomTransform.NetTranslation().x() << G4endl;
165 fout << -voxelHalfWidthY * ny + thePhantomTransform.NetTranslation().y() << " "
166 << voxelHalfWidthY * ny + thePhantomTransform.NetTranslation().y() << G4endl;
167 fout << -voxelHalfWidthZ * nz + thePhantomTransform.NetTranslation().z() << " "
168 << voxelHalfWidthZ * nz + thePhantomTransform.NetTranslation().z() << G4endl;
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 * 2 + 1) * voxelHalfWidthX,
178 (-ny + iy * 2 + 1) * voxelHalfWidthY,
179 (-nz + iz * 2 + 1) * voxelHalfWidthZ);
180 theTransform.ApplyPointTransform(voxelCentre);
181 G4bool bVoxelIsInside = true;
182 for (G4int ivx = -1; ivx <= 1; ivx += 2) {
183 for (G4int ivy = -1; ivy <= 1; ivy += 2) {
184 for (G4int ivz = -1; ivz <= 1; ivz += 2) {
185 G4ThreeVector voxelPoint = voxelCentre + ivx * voxelHalfWidthX * axisX
186 + ivy * voxelHalfWidthY * axisY
187 + ivz * voxelHalfWidthZ * axisZ;
188 if (fSolid->Inside(voxelPoint) == kOutside) {
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::SetNewValue", "", FatalException,
204 "Volume cannot intersect phantom in discontiguous voxels, "
205 "please use other voxel");
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 << G4endl;
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 * iz;
228 if (fVoxelIsInside[copyNo]) {
229 fout << thePhantomParam->GetMaterialIndex(copyNo) << " ";
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 * iz;
243 if (fVoxelIsInside[copyNo]) {
244 fout << thePhantomParam->GetMaterial(copyNo)->GetDensity() / g * cm3 << " ";
245 b1xFound = true;
246 }
247 }
248 if (b1xFound) fout << G4endl;
249 }
250 }
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 void DicomIntersectVolume::BuildUserSolid(std::vector<G4String> params)
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(tgrSolid);
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 void DicomIntersectVolume::BuildG4Solid(std::vector<G4String> params)
268 {
269 fSolid = GetLogicalVolumes(params[0], 1, 1)[0]->GetSolid();
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273 G4PhantomParameterisation* DicomIntersectVolume::GetPhantomParam(G4bool bMustExist)
274 {
275 G4PhantomParameterisation* paramreg = 0;
276
277 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
278 for (auto cite = pvs->cbegin(); cite != pvs->cend(); ++cite) {
279 if (IsPhantomVolume(*cite)) {
280 const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
281 G4VPVParameterisation* param = pvparam->GetParameterisation();
282 // if( static_cast<const G4PhantomParameterisation*>(param) ){
283 // if( static_cast<const G4PhantomParameterisation*>(param) ){
284 // G4cout << "G4PhantomParameterisation volume found "
285 // << (*cite)->GetName() << G4endl;
286 paramreg = static_cast<G4PhantomParameterisation*>(param);
287 }
288 }
289
290 if (!paramreg && bMustExist)
291 G4Exception("DicomIntersectVolume::GetPhantomParam", "", FatalErrorInArgument,
292 " No G4PhantomParameterisation found ");
293
294 return paramreg;
295 }
296
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
298 std::vector<G4VPhysicalVolume*> DicomIntersectVolume::GetPhysicalVolumes(const G4String& name,
299 G4bool exists, G4int nVols)
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(":", ial - 1);
307 if (ial2 != std::string::npos) {
308 G4Exception("DicomIntersectVolume::GetPhysicalVolumes", "", FatalErrorInArgument,
309 G4String("Name corresponds to a touchable " + name).c_str());
310 }
311 else {
312 volname = name.substr(0, ial);
313 volcopy = G4UIcommand::ConvertToInt(name.substr(ial + 1, name.length()).c_str());
314 }
315 }
316 else {
317 volname = name;
318 volcopy = -1;
319 }
320
321 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
322 for (auto citepv = pvs->cbegin(); citepv != pvs->cend(); ++citepv) {
323 if (volname == (*citepv)->GetName() && ((*citepv)->GetCopyNo() == volcopy || -1 == volcopy)) {
324 vvolu.push_back(*citepv);
325 }
326 }
327
328 //----- Check that at least one volume was found
329 if (vvolu.size() == 0) {
330 if (exists) {
331 G4Exception(" DicomIntersectVolume::GetPhysicalVolumes", "", FatalErrorInArgument,
332 G4String("No physical volume found with name " + name).c_str());
333 }
334 else {
335 G4cerr << "!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: "
336 << " no physical volume found with name " << name << G4endl;
337 }
338 }
339
340 if (nVols != -1 && G4int(vvolu.size()) != nVols) {
341 G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
342 "Wrong number of physical volumes found", FatalErrorInArgument,
343 ("Number of physical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size()))
344 + ", requesting " + G4UIcommand::ConvertToString(nVols))
345 .c_str());
346 }
347
348 return vvolu;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352 G4bool DicomIntersectVolume::IsPhantomVolume(G4VPhysicalVolume* pv)
353 {
354 EAxis axis;
355 G4int nReplicas;
356 G4double width, offset;
357 G4bool consuming;
358 pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
359 EVolume type = (consuming) ? kReplica : kParameterised;
360 if (type == kParameterised && pv->GetRegularStructureId() == 1) {
361 return TRUE;
362 }
363 else {
364 return FALSE;
365 }
366 }
367
368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
369 std::vector<G4LogicalVolume*> DicomIntersectVolume::GetLogicalVolumes(const G4String& name,
370 G4bool exists, G4int nVols)
371 {
372 // G4cout << "GetLogicalVolumes " << name << " " << exists << G4endl;
373 std::vector<G4LogicalVolume*> vvolu;
374 G4int ial = G4int(name.rfind(":"));
375 if (ial != -1) {
376 G4Exception("DicomIntersectVolume::GetLogicalVolumes", "", FatalErrorInArgument,
377 G4String("Name corresponds to a touchable or physcal volume" + name).c_str());
378 }
379
380 G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
381 for (auto citelv = lvs->cbegin(); citelv != lvs->cend(); ++citelv) {
382 if (name == (*citelv)->GetName()) {
383 vvolu.push_back(*citelv);
384 }
385 }
386
387 //----- Check that at least one volume was found
388 if (vvolu.size() == 0) {
389 if (exists) {
390 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "", FatalErrorInArgument,
391 ("no logical volume found with name " + name).c_str());
392 }
393 else {
394 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "", JustWarning,
395 ("no logical volume found with name " + name).c_str());
396 }
397 }
398
399 if (nVols != -1 && G4int(vvolu.size()) != nVols) {
400 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "Wrong number of logical volumes found",
401 FatalErrorInArgument,
402 ("Number of logical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size()))
403 + ", requesting " + G4UIcommand::ConvertToString(nVols))
404 .c_str());
405 }
406
407 return vvolu;
408 }
409
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411 std::vector<G4String> DicomIntersectVolume::GetWordsInString(const G4String& stemp)
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:GetWordsFromString", "", FatalException,
428 ("There cannot be two quotes together " + stemp).c_str());
429 }
430 if (nQuotes % 2 == 1) {
431 // close word
432 wordlist.push_back(stemp.substr(istart, ii - istart));
433 // G4cout << "GetWordsInString new word0 "
434 //<< wordlist[wordlist.size()-1] << " istart " << istart
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 << "GetWordsInString blank nQuotes "
446 //<< nQuotes << " lastIsBlank " << lastIsBlank << G4endl;
447 if (nQuotes % 2 == 0) {
448 if (!lastIsBlank && !lastIsQuote) {
449 wordlist.push_back(stemp.substr(istart, ii - istart));
450 // G4cout << "GetWordsInString new word1 "
451 //<< wordlist[wordlist.size()-1] << " lastIsBlank "
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, ii - istart + 1));
463 // G4cout << "GetWordsInString new word2 "
464 //<< wordlist[wordlist.size()-1] << " istart " << istart << G4endl;
465 }
466 lastIsQuote = false;
467 lastIsBlank = false;
468 }
469 }
470 if (nQuotes % 2 == 1) {
471 G4Exception("GmGenUtils:GetWordsFromString", "", FatalException,
472 ("unbalanced quotes in line " + stemp).c_str());
473 }
474
475 return wordlist;
476 }
477