Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
DicomIntersectVolume Class Reference

Manages intersections of DICOM files with volumes. More...

#include <DicomIntersectVolume.hh>

Inheritance diagram for DicomIntersectVolume:
G4UImessenger

Public Member Functions

 DicomIntersectVolume ()
 
 ~DicomIntersectVolume ()
 
virtual void SetNewValue (G4UIcommand *command, G4String newValues)
 
- Public Member Functions inherited from G4UImessenger
 G4UImessenger ()
 
 G4UImessenger (const G4String &path, const G4String &dsc, G4bool commandsToBeBroadcasted=true)
 
virtual ~G4UImessenger ()
 
virtual G4String GetCurrentValue (G4UIcommand *command)
 
G4bool operator== (const G4UImessenger &messenger) const
 

Additional Inherited Members

- Protected Member Functions inherited from G4UImessenger
G4String ItoS (G4int i)
 
G4String DtoS (G4double a)
 
G4String BtoS (G4bool b)
 
G4int StoI (G4String s)
 
G4double StoD (G4String s)
 
G4bool StoB (G4String s)
 
void AddUIcommand (G4UIcommand *newCommand)
 
void CreateDirectory (const G4String &path, const G4String &dsc, G4bool commandsToBeBroadcasted=true)
 
template<typename T >
T * CreateCommand (const G4String &cname, const G4String &dsc)
 
- Protected Attributes inherited from G4UImessenger
G4UIdirectory * baseDir
 
G4String baseDirName
 

Detailed Description

Manages intersections of DICOM files with volumes.

Definition at line 50 of file DicomIntersectVolume.hh.

Constructor & Destructor Documentation

DicomIntersectVolume::DicomIntersectVolume ( )

Definition at line 49 of file DicomIntersectVolume.cc.

References G4UIcommand::AvailableForStates(), G4State_Idle, G4UIcommand::SetGuidance(), and G4UIcmdWithAString::SetParameterName().

50  : G4UImessenger(),
51  fG4VolumeCmd(0),
52  fSolid(0),
53  fPhantomMinusCorner(),
54  fVoxelIsInside(0)
55 {
56  fUserVolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
57  fUserVolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume and outputs the "
58  "voxels that are totally inside the intersection as a new phantom "
59  "file. It must have the parameters: POS_X POS_Y POS_Z "
60  "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
61  fUserVolumeCmd->SetParameterName("choice",true);
62  fUserVolumeCmd->AvailableForStates(G4State_Idle);
63 
64  fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
65  fG4VolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume and outputs the "
66  "voxels that are totally inside the intersection as a new phantom "
67  "file. It must have the parameters: POS_X POS_Y POS_Z "
68  "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
69  fG4VolumeCmd->SetParameterName("choice",true);
70  fG4VolumeCmd->AvailableForStates(G4State_Idle);
71 }
G4UIcmdWithAString::SetParameterName
void SetParameterName(const char *theName, G4bool omittable, G4bool currentAsDefault=false)
Definition: G4UIcmdWithAString.cc:42
G4UIcommand::SetGuidance
void SetGuidance(const char *aGuidance)
Definition: G4UIcommand.hh:161
G4UIcommand::AvailableForStates
void AvailableForStates(G4ApplicationState s1)
Definition: G4UIcommand.cc:225
G4State_Idle
Definition: G4ApplicationState.hh:85
DicomIntersectVolume::~DicomIntersectVolume ( )

Definition at line 74 of file DicomIntersectVolume.cc.

75 {
76  if (fUserVolumeCmd) delete fUserVolumeCmd;
77  if (fG4VolumeCmd) delete fG4VolumeCmd;
78 }

Member Function Documentation

void DicomIntersectVolume::SetNewValue ( G4UIcommand *  command,
G4String  newValues 
)
virtual

Reimplemented from G4UImessenger.

Definition at line 81 of file DicomIntersectVolume.cc.

References G4AffineTransform::ApplyAxisTransform(), G4AffineTransform::ApplyPointTransform(), python.hepunit::cm3, G4UIcommand::ConvertToDouble(), G4UIcommand::ConvertToString(), FatalErrorInArgument, FatalException, g(), G4endl, G4Exception(), G4UIcommand::GetCommandName(), G4UIcommand::GetCommandPath(), G4Material::GetDensity(), G4VPhysicalVolume::GetFrameRotation(), G4VPhysicalVolume::GetFrameTranslation(), G4PhantomParameterisation::GetMaterial(), G4PhantomParameterisation::GetMaterialIndex(), G4PhantomParameterisation::GetMaterials(), G4PhantomParameterisation::GetNoVoxelX(), G4PhantomParameterisation::GetNoVoxelY(), G4PhantomParameterisation::GetNoVoxelZ(), G4PhantomParameterisation::GetVoxelHalfX(), G4PhantomParameterisation::GetVoxelHalfY(), G4PhantomParameterisation::GetVoxelHalfZ(), G4VSolid::Inside(), G4AffineTransform::Invert(), iz, kOutside, CLHEP::HepRotation::rotateX(), and CLHEP::HepRotation::rotateY().

83 {
84  G4AffineTransform theVolumeTransform;
85 
86  if (command == fUserVolumeCmd) {
87 
88  std::vector<G4String> params = GetWordsInString( newValues );
89  if( params.size() < 8 ) {
90  G4Exception("DicomIntersectVolume::SetNewValue",
91  " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
92  "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
93  FatalErrorInArgument,
94  G4String("Number of parameters given = " +
95  G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
96 
97  }
98 
99  //----- Build G4VSolid
100  BuildUserSolid(params);
101 
102  //----- Calculate volume inverse 3D transform
103  G4ThreeVector pos = G4ThreeVector( G4UIcommand::ConvertToDouble(params[0]),
104  G4UIcommand::ConvertToDouble(params[1]),
105  G4UIcommand::ConvertToDouble(params[2]) );
106  G4RotationMatrix* rotmat = new G4RotationMatrix;
107  std::vector<double> angles;
108  rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
109  rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
110  rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
111  theVolumeTransform = G4AffineTransform( rotmat, pos ).Invert();
112 
113  } else if (command == fG4VolumeCmd) {
114  std::vector<G4String> params = GetWordsInString( newValues );
115  if( params.size() !=1 ) G4Exception("DicomIntersectVolume::SetNewValue",
116  "",
117  FatalErrorInArgument,
118  G4String("Command: "+ command->GetCommandPath() +
119  "/" + command->GetCommandName() +
120  " " + newValues +
121  " needs 1 argument: VOLUME_NAME").c_str());
122 
123  //----- Build G4VSolid
124  BuildG4Solid(params);
125 
126  //----- Calculate volume inverse 3D transform
127  G4VPhysicalVolume* pv = GetPhysicalVolumes( params[0], 1, 1)[0];
128 
129  theVolumeTransform = G4AffineTransform( pv->GetFrameRotation(), pv->GetFrameTranslation() );
130  }
131 
132  //----- Calculate relative phantom - volume 3D transform
133  G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
134 
135  G4RotationMatrix* rotph = new G4RotationMatrix();
136  // assumes the phantom mother is not rotated !!!
137  G4AffineTransform thePhantomTransform( rotph, G4ThreeVector() );
138  // assumes the phantom mother is not translated !!!
139 
140  G4AffineTransform theTransform = theVolumeTransform*thePhantomTransform;
141 
142  G4ThreeVector axisX( 1., 0., 0. );
143  G4ThreeVector axisY( 0., 1., 0. );
144  G4ThreeVector axisZ( 0., 0., 1. );
145  theTransform.ApplyAxisTransform(axisX);
146  theTransform.ApplyAxisTransform(axisY);
147  theTransform.ApplyAxisTransform(axisZ);
148 
149  //----- Write phantom header
150  G4String thePhantomFileName = "phantom.g4pdcm";
151  fout.open(thePhantomFileName);
152  std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
153  fout << materials.size() << G4endl;
154  for( unsigned int ii = 0; ii < materials.size(); ii++ ) {
155  fout << ii << " " << materials[ii]->GetName() << G4endl;
156  }
157 
158  //----- Loop to pantom voxels
159  G4int nx = thePhantomParam->GetNoVoxelX();
160  G4int ny = thePhantomParam->GetNoVoxelY();
161  G4int nz = thePhantomParam->GetNoVoxelZ();
162  G4int nxy = nx*ny;
163  fVoxelIsInside = new G4bool[nx*ny*nz];
164  G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
165  G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
166  G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
167 
168  //----- Write phantom dimensions and limits
169  fout << nx << " " << ny << " " << nz << G4endl;
170  fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << " "
171  << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
172  fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << " "
173  << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
174  fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << " "
175  << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;
176 
177  for( G4int iz = 0; iz < nz; iz++ ){
178  for( G4int iy = 0; iy < ny; iy++) {
179 
180  G4bool bPrevVoxelInside = true;
181  G4bool b1VoxelFoundInside = false;
182  G4int firstVoxel = -1;
183  G4int lastVoxel = -1;
184  for(G4int ix = 0; ix < nx; ix++ ){
185  G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX,
186  (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
187  theTransform.ApplyPointTransform(voxelCentre);
188  G4bool bVoxelIsInside = true;
189  for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
190  for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
191  for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
192  G4ThreeVector voxelPoint = voxelCentre + ivx*voxelHalfWidthX*axisX +
193  ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
194  if( fSolid->Inside( voxelPoint ) == kOutside ) {
195  bVoxelIsInside = false;
196  break;
197  } else {
198  }
199  }
200  if( !bVoxelIsInside ) break;
201  }
202  if( !bVoxelIsInside ) break;
203  }
204 
205  G4int copyNo = ix + nx*iy + nxy*iz;
206  if( bVoxelIsInside ) {
207  if( !bPrevVoxelInside ) {
208  G4Exception("DicomIntersectVolume::SetNewValue",
209  "",
210  FatalException,
211  "Volume cannot intersect phantom in discontiguous voxels, "
212  "please use other voxel");
213  }
214  if( !b1VoxelFoundInside ) {
215  firstVoxel = ix;
216  b1VoxelFoundInside = true;
217  }
218  lastVoxel = ix;
219  fVoxelIsInside[copyNo] = true;
220  } else {
221  fVoxelIsInside[copyNo] = false;
222  }
223 
224  }
225  fout << firstVoxel << " " << lastVoxel << G4endl;
226  }
227  }
228 
229  //----- Now write the materials
230  for( G4int iz = 0; iz < nz; iz++ ){
231  for( G4int iy = 0; iy < ny; iy++) {
232  G4bool b1xFound = false;
233  for(G4int ix = 0; ix < nx; ix++ ){
234  size_t copyNo = ix + ny*iy + nxy*iz;
235  // fout << " iz " << iz << " i " << iy << " ix " << ix << G4endl;
236  if( fVoxelIsInside[copyNo] ) {
237  fout << thePhantomParam->GetMaterialIndex(copyNo)<< " ";
238  b1xFound = true;
239  }
240  }
241  if(b1xFound ) fout << G4endl;
242  }
243  }
244 
245  // Now write densities
246  for( G4int iz = 0; iz < nz; iz++ ){
247  for( G4int iy = 0; iy < ny; iy++) {
248  G4bool b1xFound = false;
249  for(G4int ix = 0; ix < nx; ix++ ){
250  size_t copyNo = ix + ny*iy + nxy*iz;
251  if( fVoxelIsInside[copyNo] ) {
252  fout << thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3 << " ";
253  b1xFound = true;
254  }
255  }
256  if(b1xFound ) fout << G4endl;
257  }
258  }
259 
260 }
G4VPhysicalVolume::GetFrameTranslation
G4ThreeVector GetFrameTranslation() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::HepRotation::rotateX
HepRotation & rotateX(double delta)
Definition: Rotation.cc:66
G4RotationMatrix
CLHEP::HepRotation G4RotationMatrix
Definition: G4RotationMatrix.hh:43
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:178
G4UIcommand::ConvertToString
static G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:357
CLHEP::HepRotation::rotateY
HepRotation & rotateY(double delta)
Definition: Rotation.cc:79
G4AffineTransform::ApplyAxisTransform
void ApplyAxisTransform(G4ThreeVector &axis) const
G4PhantomParameterisation::GetMaterialIndex
size_t GetMaterialIndex(size_t nx, size_t ny, size_t nz) const
Definition: G4PhantomParameterisation.cc:151
G4int
int G4int
Definition: G4Types.hh:78
G4PhantomParameterisation::GetNoVoxelZ
size_t GetNoVoxelZ() const
G4VPhysicalVolume::GetFrameRotation
const G4RotationMatrix * GetFrameRotation() const
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4AffineTransform::Invert
G4AffineTransform & Invert()
G4VSolid::Inside
virtual EInside Inside(const G4ThreeVector &p) const =0
G4bool
bool G4bool
Definition: G4Types.hh:79
iz
G4double iz
Definition: TRTMaterials.hh:39
G4PhantomParameterisation::GetMaterials
std::vector< G4Material * > GetMaterials() const
G4UIcommand::ConvertToDouble
static G4double ConvertToDouble(const char *st)
Definition: G4UIcommand.cc:429
G4UIcommand::GetCommandPath
const G4String & GetCommandPath() const
Definition: G4UIcommand.hh:139
G4PhantomParameterisation::GetVoxelHalfY
G4double GetVoxelHalfY() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhantomParameterisation::GetMaterial
G4Material * GetMaterial(size_t nx, size_t ny, size_t nz) const
Definition: G4PhantomParameterisation.cc:159
G4UIcommand::GetCommandName
const G4String & GetCommandName() const
Definition: G4UIcommand.hh:141
G4PhantomParameterisation::GetNoVoxelY
size_t GetNoVoxelY() const
FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:61
G4PhantomParameterisation::GetVoxelHalfX
G4double GetVoxelHalfX() const
FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:61
G4PhantomParameterisation::GetVoxelHalfZ
G4double GetVoxelHalfZ() const
G4PhantomParameterisation::GetNoVoxelX
size_t GetNoVoxelX() const
G4double
double G4double
Definition: G4Types.hh:76
G4AffineTransform::ApplyPointTransform
void ApplyPointTransform(G4ThreeVector &vec) const
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:01 for Geant4.10 by   doxygen 1.8.7