DicomIntersectVolume Class Reference

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

#include <DicomIntersectVolume.hh>

Inheritance diagram for DicomIntersectVolume:

Collaboration diagram for DicomIntersectVolume:

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
G4bool	CommandsShouldBeInMaster () 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
G4bool	commandsShouldBeInMaster

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.

 : G4UImessenger(),
   fG4VolumeCmd(0),
   fSolid(0),
   fPhantomMinusCorner(),
   fVoxelIsInside(0)
{
  fUserVolumeCmd= new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
  fUserVolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume "
     "and outputs the voxels that are totally inside the intersection as"
     " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
     "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
  fUserVolumeCmd->SetParameterName("choice",true);
  fUserVolumeCmd->AvailableForStates(G4State_Idle);

  fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
  fG4VolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume"
     " and outputs the voxels that are totally inside the intersection as "
     " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
     "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
  fG4VolumeCmd->SetParameterName("choice",true);
  fG4VolumeCmd->AvailableForStates(G4State_Idle);
}

Here is the call graph for this function:

DicomIntersectVolume::~DicomIntersectVolume

(

)

Definition at line 74 of file DicomIntersectVolume.cc.

{
 if (fUserVolumeCmd) delete fUserVolumeCmd;
 if (fG4VolumeCmd) delete fG4VolumeCmd;
}

Member Function Documentation

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

virtual

Reimplemented from G4UImessenger.

Definition at line 81 of file DicomIntersectVolume.cc.

{ 
  G4AffineTransform theVolumeTransform;

  if (command == fUserVolumeCmd) {

    std::vector<G4String> params = GetWordsInString( newValues );
    if( params.size() < 8 ) {
      G4Exception("DicomIntersectVolume::SetNewValue",
          " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
          "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
                  FatalErrorInArgument,
                  G4String("Number of parameters given = " + 
               G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
      
    }

    //----- Build G4VSolid 
    BuildUserSolid(params);

    //----- Calculate volume inverse 3D transform
    G4ThreeVector pos = G4ThreeVector( G4UIcommand::ConvertToDouble(params[0]), 
                       G4UIcommand::ConvertToDouble(params[1]),
                       G4UIcommand::ConvertToDouble(params[2]) );
    G4RotationMatrix* rotmat = new G4RotationMatrix;
    std::vector<double> angles; 
    rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
    rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
    rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
    theVolumeTransform = G4AffineTransform( rotmat, pos ).Invert();

  } else if (command == fG4VolumeCmd) {
    std::vector<G4String> params = GetWordsInString( newValues );
    if( params.size() !=1 ) 
      G4Exception("DicomIntersectVolume::SetNewValue",
                  "",
                  FatalErrorInArgument,
                  G4String("Command: "+ command->GetCommandPath() + 
                           "/" + command->GetCommandName() + 
                           " " + newValues + 
                           "  needs 1 argument: VOLUME_NAME").c_str()); 

    //----- Build G4VSolid 
    BuildG4Solid(params);

    //----- Calculate volume inverse 3D transform
    G4VPhysicalVolume* pv = GetPhysicalVolumes( params[0], 1, 1)[0];
    
    theVolumeTransform = 
      G4AffineTransform( pv->GetFrameRotation(), pv->GetFrameTranslation() );
  }

  //----- Calculate relative phantom - volume 3D transform 
  G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
  
  G4RotationMatrix* rotph = new G4RotationMatrix(); 
  // assumes the phantom mother is not rotated !!! 
  G4AffineTransform thePhantomTransform( rotph, G4ThreeVector() ); 
  // assumes the phantom mother is not translated !!! 

  G4AffineTransform theTransform = theVolumeTransform*thePhantomTransform;

  G4ThreeVector axisX( 1., 0., 0. );
  G4ThreeVector axisY( 0., 1., 0. );
  G4ThreeVector axisZ( 0., 0., 1. );
  theTransform.ApplyAxisTransform(axisX);
  theTransform.ApplyAxisTransform(axisY);
  theTransform.ApplyAxisTransform(axisZ);
  
  //----- Write phantom header
  G4String thePhantomFileName = "phantom.g4pdcm";
  fout.open(thePhantomFileName);
  std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
  fout << materials.size() << G4endl;
  for( unsigned int ii = 0; ii < materials.size(); ii++ ) {
    fout << ii << " " << materials[ii]->GetName() << G4endl;
  }

  //----- Loop to pantom voxels
  G4int nx = thePhantomParam->GetNoVoxelX();
  G4int ny = thePhantomParam->GetNoVoxelY();
  G4int nz = thePhantomParam->GetNoVoxelZ();
  G4int nxy = nx*ny;
  fVoxelIsInside = new G4bool[nx*ny*nz];
  G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
  G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
  G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();

  //----- Write phantom dimensions and limits
  fout << nx << " " << ny << " " << nz << G4endl;
  fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << " "  
       << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
  fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << " " 
       << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
  fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << " " 
       << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;

  for( G4int iz = 0; iz < nz; iz++ ){
    for( G4int iy = 0; iy < ny; iy++) {

      G4bool bPrevVoxelInside = true;
      G4bool b1VoxelFoundInside = false;
      G4int firstVoxel = -1;
      G4int lastVoxel = -1;
      for(G4int ix = 0; ix < nx; ix++ ){
        G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX, 
                   (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
        theTransform.ApplyPointTransform(voxelCentre);
        G4bool bVoxelIsInside = true;
        for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
          for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
            for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
              G4ThreeVector voxelPoint = voxelCentre 
                + ivx*voxelHalfWidthX*axisX +
                ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
              if( fSolid->Inside( voxelPoint ) == kOutside ) {
                bVoxelIsInside = false;
                break;
              } else {
              }
            }
            if( !bVoxelIsInside ) break;
          }
          if( !bVoxelIsInside ) break;
        }

        G4int copyNo = ix + nx*iy + nxy*iz;
        if( bVoxelIsInside ) {
          if( !bPrevVoxelInside ) {
            G4Exception("DicomIntersectVolume::SetNewValue",
                        "",
                        FatalException,
                 "Volume cannot intersect phantom in discontiguous voxels, "
                 "please use other voxel");
          }
          if( !b1VoxelFoundInside ) {
            firstVoxel = ix;
            b1VoxelFoundInside = true;
          }
          lastVoxel = ix;
          fVoxelIsInside[copyNo] = true;
        } else {
          fVoxelIsInside[copyNo] = false;
        }
           
      }
      fout << firstVoxel << " " << lastVoxel << G4endl;
    }
  }

  //-----  Now write the materials
  for( G4int iz = 0; iz < nz; iz++ ){
    for( G4int iy = 0; iy < ny; iy++) {
      G4bool b1xFound = false; 
      for(G4int ix = 0; ix < nx; ix++ ){
        size_t copyNo = ix + ny*iy + nxy*iz;
        //        fout << " iz " << iz << " i " << iy << " ix " << ix << G4endl;
        if( fVoxelIsInside[copyNo] ) {
          fout << thePhantomParam->GetMaterialIndex(copyNo)<< " ";
          b1xFound = true;
        }
      }
      if(b1xFound ) fout << G4endl;
    }
  }

  // Now write densities
  for( G4int iz = 0; iz < nz; iz++ ){
    for( G4int iy = 0; iy < ny; iy++) {
      G4bool b1xFound = false; 
      for(G4int ix = 0; ix < nx; ix++ ){
        size_t copyNo = ix + ny*iy + nxy*iz;
        if( fVoxelIsInside[copyNo] ) {
          fout <<thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3<< " ";
          b1xFound = true;
        }
      }
      if(b1xFound ) fout << G4endl;
    }
  }

} 

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• source/geant4.10.03.p03/examples/extended/medical/DICOM/include/DicomIntersectVolume.hh
• source/geant4.10.03.p03/examples/extended/medical/DICOM/src/DicomIntersectVolume.cc