ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomIntersectVolume.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DicomIntersectVolume.cc
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 //
29 //
30 
31 #include "DicomIntersectVolume.hh"
32 
33 #include "G4UIcmdWithAString.hh"
34 #include "G4LogicalVolume.hh"
35 #include "G4LogicalVolumeStore.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4PhysicalVolumeStore.hh"
38 #include "G4VSolid.hh"
39 #include "G4tgrSolid.hh"
40 #include "G4tgbVolume.hh"
41 #include "G4Material.hh"
43 #include "G4PVParameterised.hh"
44 #include "G4UIcommand.hh"
45 #include "G4SystemOfUnits.hh"
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49  : G4UImessenger(),
50  fG4VolumeCmd(0),
51  fSolid(0),
52  fPhantomMinusCorner(),
53  fVoxelIsInside(0)
54 {
55  fUserVolumeCmd= new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
56  fUserVolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume "
57  "and outputs the voxels that are totally inside the intersection as"
58  " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
59  "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
60  fUserVolumeCmd->SetParameterName("choice",true);
62 
63  fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
64  fG4VolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume"
65  " and outputs the voxels that are totally inside the intersection as "
66  " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
67  "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
68  fG4VolumeCmd->SetParameterName("choice",true);
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 {
75  if (fUserVolumeCmd) delete fUserVolumeCmd;
76  if (fG4VolumeCmd) delete fG4VolumeCmd;
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81  G4String newValues)
82 {
83  G4AffineTransform theVolumeTransform;
84 
85  if (command == fUserVolumeCmd)
86  {
87  std::vector<G4String> params = GetWordsInString( newValues );
88  if( params.size() < 8 ) {
89  G4Exception("DicomIntersectVolume::SetNewValue",
90  " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
91  "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
93  G4String("Number of parameters given = " +
94  G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
95 
96  }
97 
98  //----- Build G4VSolid
99  BuildUserSolid(params);
100 
101  //----- Calculate volume inverse 3D transform
103  G4UIcommand::ConvertToDouble(params[1]),
104  G4UIcommand::ConvertToDouble(params[2]) );
105  G4RotationMatrix* rotmat = new G4RotationMatrix;
106  std::vector<G4double> angles;
107  rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
108  rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
109  rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
110  theVolumeTransform = G4AffineTransform( rotmat, pos ).Invert();
111 
112  }
113  else if (command == fG4VolumeCmd)
114  {
115  std::vector<G4String> params = GetWordsInString( newValues );
116  if( params.size() !=1 )
117  G4Exception("DicomIntersectVolume::SetNewValue",
118  "",
120  G4String("Command: "+ command->GetCommandPath() +
121  "/" + command->GetCommandName() +
122  " " + newValues +
123  " needs 1 argument: VOLUME_NAME").c_str());
124 
125  //----- Build G4VSolid
126  BuildG4Solid(params);
127 
128  //----- Calculate volume inverse 3D transform
129  G4VPhysicalVolume* pv = GetPhysicalVolumes( params[0], 1, 1)[0];
130 
131  theVolumeTransform =
133  }
134 
135  //----- Calculate relative phantom - volume 3D transform
136  G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
137 
138  G4RotationMatrix* rotph = new G4RotationMatrix();
139  // assumes the phantom mother is not rotated !!!
140  G4AffineTransform thePhantomTransform( rotph, G4ThreeVector() );
141  // assumes the phantom mother is not translated !!!
142 
143  G4AffineTransform theTransform = theVolumeTransform*thePhantomTransform;
144 
145  G4ThreeVector axisX( 1., 0., 0. );
146  G4ThreeVector axisY( 0., 1., 0. );
147  G4ThreeVector axisZ( 0., 0., 1. );
148  theTransform.ApplyAxisTransform(axisX);
149  theTransform.ApplyAxisTransform(axisY);
150  theTransform.ApplyAxisTransform(axisZ);
151 
152  //----- Write phantom header
153  G4String thePhantomFileName = "phantom.g4pdcm";
154  fout.open(thePhantomFileName);
155  std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
156  fout << materials.size() << G4endl;
157  for( unsigned int ii = 0; ii < materials.size(); ++ii )
158  {
159  fout << ii << " " << materials[ii]->GetName() << G4endl;
160  }
161 
162  //----- Loop to pantom voxels
163  G4int nx = G4int(thePhantomParam->GetNoVoxelX());
164  G4int ny = G4int(thePhantomParam->GetNoVoxelY());
165  G4int nz = G4int(thePhantomParam->GetNoVoxelZ());
166  G4int nxy = nx*ny;
167  fVoxelIsInside = new G4bool[nx*ny*nz];
168  G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
169  G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
170  G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
171 
172  //----- Write phantom dimensions and limits
173  fout << nx << " " << ny << " " << nz << G4endl;
174  fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << " "
175  << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
176  fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << " "
177  << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
178  fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << " "
179  << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;
180 
181  for( G4int iz = 0; iz < nz; ++iz ){
182  for( G4int iy = 0; iy < ny; ++iy) {
183 
184  G4bool bPrevVoxelInside = true;
185  G4bool b1VoxelFoundInside = false;
186  G4int firstVoxel = -1;
187  G4int lastVoxel = -1;
188  for(G4int ix = 0; ix < nx; ++ix ){
189  G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX,
190  (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
191  theTransform.ApplyPointTransform(voxelCentre);
192  G4bool bVoxelIsInside = true;
193  for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
194  for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
195  for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
196  G4ThreeVector voxelPoint = voxelCentre
197  + ivx*voxelHalfWidthX*axisX +
198  ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
199  if( fSolid->Inside( voxelPoint ) == kOutside ) {
200  bVoxelIsInside = false;
201  break;
202  } else {
203  }
204  }
205  if( !bVoxelIsInside ) break;
206  }
207  if( !bVoxelIsInside ) break;
208  }
209 
210  G4int copyNo = ix + nx*iy + nxy*iz;
211  if( bVoxelIsInside ) {
212  if( !bPrevVoxelInside ) {
213  G4Exception("DicomIntersectVolume::SetNewValue",
214  "",
216  "Volume cannot intersect phantom in discontiguous voxels, "
217  "please use other voxel");
218  }
219  if( !b1VoxelFoundInside ) {
220  firstVoxel = ix;
221  b1VoxelFoundInside = true;
222  }
223  lastVoxel = ix;
224  fVoxelIsInside[copyNo] = true;
225  } else {
226  fVoxelIsInside[copyNo] = false;
227  }
228 
229  }
230  fout << firstVoxel << " " << lastVoxel << G4endl;
231  }
232  }
233 
234  //----- Now write the materials
235  for( G4int iz = 0; iz < nz; ++iz ){
236  for( G4int iy = 0; iy < ny; ++iy) {
237  G4bool b1xFound = false;
238  for(G4int ix = 0; ix < nx; ++ix ){
239  size_t copyNo = ix + ny*iy + nxy*iz;
240  if( fVoxelIsInside[copyNo] ) {
241  fout << thePhantomParam->GetMaterialIndex(copyNo)<< " ";
242  b1xFound = true;
243  }
244  }
245  if(b1xFound ) fout << G4endl;
246  }
247  }
248 
249  // Now write densities
250  for( G4int iz = 0; iz < nz; ++iz ){
251  for( G4int iy = 0; iy < ny; ++iy) {
252  G4bool b1xFound = false;
253  for(G4int ix = 0; ix < nx; ++ix ){
254  size_t copyNo = ix + ny*iy + nxy*iz;
255  if( fVoxelIsInside[copyNo] ) {
256  fout <<thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3<< " ";
257  b1xFound = true;
258  }
259  }
260  if(b1xFound ) fout << G4endl;
261  }
262  }
263 
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 void DicomIntersectVolume::BuildUserSolid( std::vector<G4String> params )
268 {
269  for( G4int ii = 0; ii < 6; ++ii ) params.erase( params.begin() );
270  // take otu position and rotation angles
271  params.insert( params.begin(), ":SOLID");
272  params.insert( params.begin(), params[1] );
273  G4tgrSolid* tgrSolid = new G4tgrSolid(params);
274  G4tgbVolume* tgbVolume = new G4tgbVolume();
275  fSolid = tgbVolume->FindOrConstructG4Solid( tgrSolid );
276 
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280 void DicomIntersectVolume::BuildG4Solid( std::vector<G4String> params )
281 {
282  fSolid = GetLogicalVolumes( params[0], 1, 1)[0]->GetSolid();
283 
284 }
285 
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 {
290  G4PhantomParameterisation* paramreg = 0;
291 
293  for( auto cite = pvs->cbegin(); cite != pvs->cend(); ++cite )
294  {
295  if( IsPhantomVolume( *cite ) ) {
296  const G4PVParameterised* pvparam =
297  static_cast<const G4PVParameterised*>(*cite);
298  G4VPVParameterisation* param = pvparam->GetParameterisation();
299  // if( static_cast<const G4PhantomParameterisation*>(param) ){
300  // if( static_cast<const G4PhantomParameterisation*>(param) ){
301  // G4cout << "G4PhantomParameterisation volume found "
302  // << (*cite)->GetName() << G4endl;
303  paramreg = static_cast<G4PhantomParameterisation*>(param);
304  }
305  }
306 
307  if( !paramreg && bMustExist )
308  G4Exception("DicomIntersectVolume::GetPhantomParam",
309  "",
311  " No G4PhantomParameterisation found ");
312 
313  return paramreg;
314 
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
318 std::vector<G4VPhysicalVolume*> DicomIntersectVolume::GetPhysicalVolumes(
319  const G4String& name, G4bool exists, G4int nVols )
320 {
321  std::vector<G4VPhysicalVolume*> vvolu;
322  std::string::size_type ial = name.rfind(":");
323  G4String volname = "";
324  G4int volcopy = 0;
325  if( ial != std::string::npos )
326  {
327  std::string::size_type ial2 = name.rfind(":",ial-1);
328  if( ial2 != std::string::npos )
329  {
330  G4Exception("DicomIntersectVolume::GetPhysicalVolumes",
331  "",
333  G4String("Name corresponds to a touchable " + name).c_str());
334  }
335  else
336  {
337  volname = name.substr( 0, ial );
338  volcopy = G4UIcommand::ConvertToInt( name.substr(ial+1, name.length()).c_str() );
339  }
340  }
341  else
342  {
343  volname = name;
344  volcopy = -1;
345  }
346 
348  for( auto citepv = pvs->cbegin(); citepv != pvs->cend(); ++citepv )
349  {
350  if( volname == (*citepv)->GetName()
351  && ( (*citepv)->GetCopyNo() == volcopy || -1 == volcopy ) )
352  {
353  vvolu.push_back( *citepv );
354  }
355  }
356 
357  //----- Check that at least one volume was found
358  if( vvolu.size() == 0 ) {
359  if(exists) {
360  G4Exception(" DicomIntersectVolume::GetPhysicalVolumes",
361  "",
363  G4String("No physical volume found with name " + name).c_str());
364  } else {
365  G4cerr << "!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: "<<
366  " no physical volume found with name " << name << G4endl;
367  }
368  }
369 
370  if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
371  G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
372  "Wrong number of physical volumes found",
374  ("Number of physical volumes " +
375  G4UIcommand::ConvertToString(G4int(vvolu.size())) +
376  ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
377  }
378 
379  return vvolu;
380 }
381 
382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
384 {
385  EAxis axis;
386  G4int nReplicas;
388  G4bool consuming;
389  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
390  EVolume type = (consuming) ? kReplica : kParameterised;
391  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {
392  return TRUE;
393  } else {
394  return FALSE;
395  }
396 
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400 std::vector<G4LogicalVolume*> DicomIntersectVolume::GetLogicalVolumes(
401 const G4String& name, G4bool exists, G4int nVols )
402 {
403  // G4cout << "GetLogicalVolumes " << name << " " << exists << G4endl;
404  std::vector<G4LogicalVolume*> vvolu;
405  G4int ial = G4int(name.rfind(":"));
406  if( ial != -1 ) {
407  G4Exception("DicomIntersectVolume::GetLogicalVolumes",
408  "",
410  G4String("Name corresponds to a touchable or physcal volume"
411  + name).c_str());
412  }
413 
415  for( auto citelv = lvs->cbegin(); citelv != lvs->cend(); ++citelv )
416  {
417  if( name == (*citelv)->GetName() ) {
418  vvolu.push_back( *citelv );
419  }
420  }
421 
422  //----- Check that at least one volume was found
423  if( vvolu.size() == 0 ) {
424  if(exists) {
425  G4Exception("DicomIntersectVolume::GetLogicalVolumes:","",
426  FatalErrorInArgument,("no logical volume found with name "
427  + name).c_str());
428  } else {
429  G4Exception("DicomIntersectVolume::GetLogicalVolumes:","",
430  JustWarning,("no logical volume found with name " + name).c_str());
431  }
432  }
433 
434  if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
435  G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
436  "Wrong number of logical volumes found",
438  ("Number of logical volumes " +
439  G4UIcommand::ConvertToString(G4int(vvolu.size())) +
440  ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
441  }
442 
443  return vvolu;
444 
445 }
446 
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
449  const G4String& stemp)
450 {
451  std::vector<G4String> wordlist;
452 
453  //---------- Read a line of file:
454  //----- Clear wordlist
455  G4int ii;
456  const char* cstr = stemp.c_str();
457  G4int siz = G4int(stemp.length());
458  G4int istart = 0;
459  G4int nQuotes = 0;
460  G4bool lastIsBlank = false;
461  G4bool lastIsQuote = false;
462  for( ii = 0; ii < siz; ++ii )
463  {
464  if(cstr[ii] == '\"' ){
465  if( lastIsQuote ){
466  G4Exception("GmGenUtils:GetWordsFromString","",FatalException,
467  ("There cannot be two quotes together " + stemp).c_str() );
468  }
469  if( nQuotes%2 == 1 ){
470  //close word
471  wordlist.push_back( stemp.substr(istart,ii-istart) );
472  // G4cout << "GetWordsInString new word0 "
473  //<< wordlist[wordlist.size()-1] << " istart " << istart
474  //<< " ii " << ii << G4endl;
475  } else {
476  istart = ii+1;
477  }
478  ++nQuotes;
479  lastIsQuote = true;
480  lastIsBlank = false;
481  } else if(cstr[ii] == ' ' ){
482  // G4cout << "GetWordsInString blank nQuotes "
483  //<< nQuotes << " lastIsBlank " << lastIsBlank << G4endl;
484  if( nQuotes%2 == 0 ){
485  if( !lastIsBlank && !lastIsQuote ) {
486  wordlist.push_back( stemp.substr(istart,ii-istart) );
487  // G4cout << "GetWordsInString new word1 "
488  //<< wordlist[wordlist.size()-1] << " lastIsBlank "
489  //<< lastIsBlank << G4endl;
490  }
491 
492  istart = ii+1;
493  lastIsQuote = false;
494  lastIsBlank = true;
495  }
496  } else {
497  if( ii == siz-1 ) {
498  wordlist.push_back( stemp.substr(istart,ii-istart+1) );
499  // G4cout << "GetWordsInString new word2 "
500  //<< wordlist[wordlist.size()-1] << " istart " << istart << G4endl;
501  }
502  lastIsQuote = false;
503  lastIsBlank = false;
504  }
505  }
506  if( nQuotes%2 == 1 ) {
507  G4Exception("GmGenUtils:GetWordsFromString","",FatalException,
508  ("unbalanced quotes in line " + stemp).c_str() );
509  }
510 
511  return wordlist;
512 }