Construct a DICOM Geant4 geometry produced from the intersetion of a DICOM file and a volume. More...

#include <Doxymodules_medical.h>

Inheritance diagram for DicomPartialDetectorConstruction:

Public Member Functions
	DicomPartialDetectorConstruction ()

	~DicomPartialDetectorConstruction ()

virtual G4VPhysicalVolume *	Construct ()

Public Member Functions inherited from DicomDetectorConstruction
	DicomDetectorConstruction ()

	~DicomDetectorConstruction ()

G4int	GetTotalVoxels () const

Private Member Functions
virtual void	ReadPhantomData ()

void	ReadPhantomDataFile (const G4String &fname)

void	ConstructPhantomContainer ()

virtual void	ConstructPhantom ()

void	ReadVoxelDensitiesPartial (std::ifstream &fin)

Private Attributes
G4PartialPhantomParameterisation *	fPartialPhantomParam

std::multimap< G4int, G4int >	fFilledIDs

std::map< G4int, std::map< G4int, G4int > >	fFilledMins

std::map< G4int, std::map< G4int, G4int > >	fFilledMaxs

G4int	fNoVoxels

G4double	fDimX

G4double	fDimY

G4double	fDimZ

G4double	fOffsetX

G4double	fOffsetY

G4double	fOffsetZ

std::vector< G4Material * >	fPhantomMaterials

Additional Inherited Members
Protected Member Functions inherited from DicomDetectorConstruction
void	InitialisationOfMaterials ()

void	ReadPhantomData ()

void	ReadPhantomDataNew ()

void	ReadVoxelDensities (std::ifstream &fin)

void	ReadPhantomDataFile (const G4String &fname)

void	MergeZSliceHeaders ()

G4Material *	BuildMaterialWithChangingDensity (const G4Material *origMate, G4float density, G4String newMateName)

void	ConstructPhantomContainer ()

void	ConstructPhantomContainerNew ()

void	SetScorer (G4LogicalVolume *voxel_logic)

virtual void	ConstructSDandField ()

Protected Attributes inherited from DicomDetectorConstruction
G4Material *	fAir

G4Box *	fWorld_solid

G4LogicalVolume *	fWorld_logic

G4VPhysicalVolume *	fWorld_phys

G4Box *	fContainer_solid

G4LogicalVolume *	fContainer_logic

G4VPhysicalVolume *	fContainer_phys

G4int	fNoFiles

std::vector< G4Material * >	fOriginalMaterials

std::vector< G4Material * >	fMaterials

size_t *	fMateIDs

std::map< G4int, G4double >	fDensityDiffs

std::vector< DicomPhantomZSliceHeader * >	fZSliceHeaders

DicomPhantomZSliceHeader *	fZSliceHeaderMerged

G4int	fNoVoxelsX

G4int	fNoVoxelsY

G4int	fNoVoxelsZ

G4double	fVoxelHalfDimX

G4double	fVoxelHalfDimY

G4double	fVoxelHalfDimZ

G4double	fMinX

G4double	fMinY

G4double	fMinZ

G4double	fMaxX

G4double	fMaxY

G4double	fMaxZ

std::map< G4int, G4Material * >	fPhantomMaterialsOriginal

DicomPhantomZSliceMerged *	fMergedSlices

std::set< G4LogicalVolume * >	fScorers

G4bool	fConstructed

Detailed Description

Construct a DICOM Geant4 geometry produced from the intersetion of a DICOM file and a volume.

Definition at line 25 of file Doxymodules_medical.h.

Constructor & Destructor Documentation

◆ DicomPartialDetectorConstruction()

DicomPartialDetectorConstruction::DicomPartialDetectorConstruction ( )

Definition at line 58 of file DicomPartialDetectorConstruction.cc.

 : DicomDetectorConstruction(),
   fPartialPhantomParam(0),
   fNoVoxels(0),
   fDimX(0),
   fDimY(0),
   fDimZ(0),
   fOffsetX(0),
   fOffsetY(0),
   fOffsetZ(0)
{
}

◆ ~DicomPartialDetectorConstruction()

DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction ( )

Definition at line 72 of file DicomPartialDetectorConstruction.cc.

73{

74}

Member Function Documentation

◆ Construct()

G4VPhysicalVolume * DicomPartialDetectorConstruction::Construct ( )

virtual

Reimplemented from DicomDetectorConstruction.

Definition at line 77 of file DicomPartialDetectorConstruction.cc.

{
  InitialisationOfMaterials();
 
  //----- Build world
  G4double worldXDimension = 1.*m;
  G4double worldYDimension = 1.*m;
  G4double worldZDimension = 1.*m;
 
  G4Box* world_box = new G4Box( "WorldSolid",
                          worldXDimension,
                          worldYDimension,
                          worldZDimension );
 
  fWorld_logic = new G4LogicalVolume( world_box,
                                    fAir,
                                    "WorldLogical",
                                    0, 0, 0 );
 
  G4VPhysicalVolume* world_phys = new G4PVPlacement( 0,
                                  G4ThreeVector(0,0,0),
                                  "World",
                                  fWorld_logic,
                                  0,
                                  false,
                                  0 );
 
  ReadPhantomData();
 
  ConstructPhantomContainer();
  ConstructPhantom();
 
  return world_phys;
}

◆ ReadPhantomData()

void DicomPartialDetectorConstruction::ReadPhantomData ( )

privatevirtual

Definition at line 113 of file DicomPartialDetectorConstruction.cc.

{
    G4String dataFile = "Data.dat";
    std::ifstream finDF(dataFile.c_str());
    G4String fname;
    if(finDF.good() != 1 ) {
        G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData",
                    "",
                    FatalErrorInArgument,
                    " Problem reading data file: Data.dat");
    }
 
    G4int compression;
    finDF >> compression; // not used here
 
    G4int numFiles;
    finDF >> numFiles;  // only 1 file supported for the moment
    if( numFiles != 1 ) {
        G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
                    "",
                    FatalErrorInArgument,
                    "More than 1 DICOM file is not supported");
    }
    for(G4int i = 0; i < numFiles; i++ ) {
        finDF >> fname;
        //--- Read one data file
        fname += ".g4pdcm";
        ReadPhantomDataFile(fname);
    }
 
    finDF.close();
}

◆ ReadPhantomDataFile()

void DicomPartialDetectorConstruction::ReadPhantomDataFile ( const G4String & fname )

private

Definition at line 147 of file DicomPartialDetectorConstruction.cc.

{
  //  G4String filename = "phantom.g4pdcm";
#ifdef G4VERBOSE
  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
         << fname << G4endl;
#endif
  std::ifstream fin(fname.c_str(), std::ios_base::in);
  if( !fin.is_open() ) {
    G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile",
                "",
                FatalException,
                G4String("File not found " + fname).c_str());
  }
  G4int nMaterials;
  fin >> nMaterials;
  G4String stemp;
  G4int nmate;
  for( G4int ii = 0; ii < nMaterials; ii++ ){
    fin >> nmate >> stemp;
    G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate "
    << ii << " = " << nmate << " mate " << stemp << G4endl;
    if( ii != nmate )
      G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
                  "",
                  FatalErrorInArgument,
                  "Material number should be in increasing order: \
                          wrong material number ");
 
  }
 
  fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z "
  << fNoVoxelsX << " " << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
  fin >> fOffsetX >> fDimX;
  fDimX = (fDimX - fOffsetX)/fNoVoxelsX;
  fin >> fOffsetY >> fDimY;
  fDimY = (fDimY - fOffsetY)/fNoVoxelsY;
  fin >> fOffsetZ >> fDimZ;
  fDimZ = (fDimZ - fOffsetZ)/fNoVoxelsZ;
  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimX "
  << fDimX << " fOffsetX " << fOffsetX << G4endl;
  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY "
  << fDimY << " fOffsetY " << fOffsetY << G4endl;
  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ "
  << fDimZ << " fOffsetZ " << fOffsetZ << G4endl;
 
  //--- Read voxels that are filled
  fNoVoxels = 0;
  //  G4bool* isFilled = new G4bool[fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ];
  //  fFilledIDs = new size_t[fNoVoxelsZ*fNoVoxelsY+1];
  //?  fFilledIDs.insert(0);
  G4int ifxmin1, ifxmax1;
  for( G4int iz = 0; iz < fNoVoxelsZ; iz++ ) {
    std::map< G4int, G4int > ifmin, ifmax;
    for( G4int iy = 0; iy < fNoVoxelsY; iy++ ) {
      fin >> ifxmin1 >> ifxmax1;
      // check coherence ...
 
      ifmin[iy] = ifxmin1;
      ifmax[iy] = ifxmax1;
      G4int ifxdiff = ifxmax1-ifxmin1+1;
      if( ifxmax1 == -1 && ifxmin1 == -1 ) ifxdiff = 0;
      fFilledIDs.insert(std::pair<G4int,G4int>(ifxdiff+fNoVoxels-1, ifxmin1));
      G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert "
             << " FilledIDs "  << ifxdiff+fNoVoxels-1 << " min " << ifxmin1
             << " N= " << fFilledIDs.size() << G4endl;
      //filledIDs[iz*fNoVoxelsY+iy+1] = ifxmax1-ifxmin1+1;
      for( G4int ix = 0; ix < fNoVoxelsX; ix++ ) {
        if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
          fNoVoxels++;
        } else {
        }
      }
    }
    fFilledMins[iz] = ifmin;
    fFilledMaxs[iz] = ifmax;
  }
 
  //--- Read material IDs
  G4int mateID1;
  fMateIDs = new size_t[fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ];
  G4int copyNo = 0;
  for( G4int iz = 0; iz < fNoVoxelsZ; iz++ ) {
    std::map< G4int, G4int > ifmin = fFilledMins[iz];
    std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
    for( G4int iy = 0; iy < fNoVoxelsY; iy++ ) {
      ifxmin1 = ifmin[iy];
      ifxmax1 = ifmax[iy];
      for( G4int ix = 0; ix < fNoVoxelsX; ix++ ) {
        if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
          fin >> mateID1;
          if( mateID1 < 0 || mateID1 >= nMaterials ) {
            G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
                        "",
                        FatalException,
            G4String("Wrong index in phantom file: It should be between 0 and "
                     + G4UIcommand::ConvertToString(nMaterials-1)
                     + ", while it is "
                     + G4UIcommand::ConvertToString(mateID1)).c_str());
          }
          fMateIDs[copyNo] = mateID1;
          copyNo++;
        }
      }
    }
  }
 
  ReadVoxelDensitiesPartial( fin );
 
  fin.close();
}

◆ ConstructPhantomContainer()

void DicomPartialDetectorConstruction::ConstructPhantomContainer ( )

private

Definition at line 407 of file DicomPartialDetectorConstruction.cc.

{
  //build a clinder that encompass the partial phantom built in the example
  // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
  //---- Extract number of voxels and voxel dimensions
 
  G4String contname= "PhantomContainer";
 
  //----- Define the volume that contains all the voxels
  G4Tubs* container_tube = new G4Tubs(contname,0., 200., 100., 0., 360*deg);
 
  fContainer_logic =
    new G4LogicalVolume( container_tube,
                         fAir,
                         contname,
                         0, 0, 0 );
 
  G4ThreeVector posCentreVoxels(0.,0.,0.);
  //G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
  G4RotationMatrix* rotm = new G4RotationMatrix;
 
  fContainer_phys =
    new G4PVPlacement(rotm,  // rotation
                      posCentreVoxels,
                      fContainer_logic,     // The logic volume
                      contname,  // Name
                      fWorld_logic,  // Mother
                      false,           // No op. bool.
                      1);              // Copy number
 
}

◆ ConstructPhantom()

void DicomPartialDetectorConstruction::ConstructPhantom ( )

privatevirtual

Implements DicomDetectorConstruction.

Definition at line 440 of file DicomPartialDetectorConstruction.cc.

{
  G4String OptimAxis = "kXAxis";
  G4bool bSkipEqualMaterials = 0;
  G4int regStructureID = 2;
 
  G4ThreeVector posCentreVoxels(0.,0.,0.);
 
  //----- Build phantom
  G4String voxelName = "phantom";
  G4Box* phantom_solid = new G4Box(voxelName, fDimX/2., fDimY/2., fDimZ/2. );
  G4LogicalVolume* phantom_logic =
    new G4LogicalVolume( phantom_solid,
                         fAir,
                         voxelName,
                         0, 0, 0 );
  G4bool pVis = 0;
  if( !pVis ) {
    G4VisAttributes* visAtt = new G4VisAttributes( FALSE );
    phantom_logic->SetVisAttributes( visAtt );
  }
 
  G4double theSmartless = 0.5;
  //  fContainer_logic->SetSmartless( theSmartless );
  phantom_logic->SetSmartless( theSmartless );
 
  fPartialPhantomParam = new G4PartialPhantomParameterisation();
 
  fPartialPhantomParam->SetMaterials( fPhantomMaterials );
  fPartialPhantomParam->SetVoxelDimensions( fDimX/2., fDimY/2., fDimZ/2. );
  fPartialPhantomParam->SetNoVoxels( fNoVoxelsX, fNoVoxelsY, fNoVoxelsZ );
  fPartialPhantomParam->SetMaterialIndices( fMateIDs );
 
  fPartialPhantomParam->SetFilledIDs(fFilledIDs);
 
  fPartialPhantomParam->SetFilledMins(fFilledMins);
 
  fPartialPhantomParam->BuildContainerWalls();
 
  //  G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
  //  G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
 
  G4PVParameterised * phantom_phys = 0;
  if( OptimAxis == "kUndefined" ) {
    phantom_phys =
      new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
                            kUndefined, fNoVoxels, fPartialPhantomParam);
  } else   if( OptimAxis == "kXAxis" ) {
    //    G4cout << " optim kX " << G4endl;
    phantom_phys =
      new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
                            kXAxis, fNoVoxels, fPartialPhantomParam);
    fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
  } else {
    G4Exception("GmReadPhantomGeometry::ConstructPhantom",
                "",
                FatalErrorInArgument,
                G4String("Wrong argument to parameter \
         GmReadPhantomGeometry:Phantom:OptimAxis: \
         Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str());
  }
 
  phantom_phys->SetRegularStructureId(regStructureID);
 
}

◆ ReadVoxelDensitiesPartial()

void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial ( std::ifstream & fin )

private

Definition at line 262 of file DicomPartialDetectorConstruction.cc.

{
    std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
    std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
    for( G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii )
    {
       densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
    }
 
    //----- Define density differences (maximum density difference to create
    // a new material)
    char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
    G4double densityDiff = -1.;
    if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
    std::map<G4int,G4double> densityDiffs;
    if( densityDiff != -1. )
    {
        for( G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii )
        {
            densityDiffs[ii] = densityDiff; //currently all materials
            //with same step
        }
    }
    else
    {
        if( fMaterials.size() == 0 ) { // do it only for first slice
            for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
            {
                fMaterials.push_back( fOriginalMaterials[ii] );
            }
        }
    }
  //  densityDiffs[0] = 0.0001; //fAir
 
  //--- Calculate the average material density for each material/density bin
  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
 
  G4double dens1;
  G4int ifxmin1, ifxmax1;
 
  //---- Read the material densities
  G4int copyNo = 0;
  for( G4int iz = 0; iz < fNoVoxelsZ; ++iz )
  {
    std::map< G4int, G4int > ifmin = fFilledMins[iz];
    std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
    for( G4int iy = 0; iy < fNoVoxelsY; ++iy )
    {
      ifxmin1 = ifmin[iy];
      ifxmax1 = ifmax[iy];
      for( G4int ix = 0; ix < fNoVoxelsX; ++ix )
      {
        if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
          fin >> dens1;
        //--- store the minimum and maximum density for each material
          //(just for printing)
          mpite = densiMinMax.find( G4int(fMateIDs[copyNo]) );
          if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
          if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
 
          //--- Get material from original list of material in file
          G4int mateID = G4int(fMateIDs[copyNo]);
          G4Material* mate = fOriginalMaterials[mateID];
          //        G4cout << copyNo << " mateID " << mateID << G4endl;
          //--- Check if density is equal to the original material density
          if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 ) continue;
 
          //--- Build material name with fOriginalMaterials name + density
          //        float densityBin = densityDiffs[mateID] *
          //        (G4int(dens1/densityDiffs[mateID])+0.5);
          G4String newMateName = mate->GetName();
 
          G4int densityBin = 0;
          //        G4cout << " densityBin " << densityBin << " "
          //        << dens1 << " " <<densityDiffs[mateID] << G4endl;
 
          if( densityDiff != -1.) {
            densityBin = (G4int(dens1/densityDiffs[mateID]));
            newMateName =
              mate->GetName()+G4UIcommand::ConvertToString(densityBin);
          }
 
          //--- Look if it is the first voxel with this material/densityBin
          std::pair<G4Material*,G4int> matdens(mate, densityBin );
 
          auto mppite = newMateDens.find( matdens );
          if( mppite != newMateDens.cend() ){
            matInfo* mi = (*mppite).second;
            mi->fSumdens += dens1;
            mi->fNvoxels++;
            fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->fId;
          //  G4cout << copyNo << " mat new again "
          //<< fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
          } else {
            matInfo* mi = new matInfo;
            mi->fSumdens = dens1;
            mi->fNvoxels = 1;
            mi->fId = G4int(newMateDens.size()+1);
            newMateDens[matdens] = mi;
            fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->fId;
            //          G4cout << copyNo << " mat new first "
            //          << fOriginalMaterials.size()-1 + mi->fId << G4endl;
          }
          copyNo++;
        //        G4cout << ix << " " << iy << " " << iz
        //  << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
                                      //        fMateIDs[copyNo] = atoi(cid)-1;
        }
      }
    }
  }
 
  //----- Build the list of phantom materials that go to Parameterisation
  //--- Add original materials
  for( auto mimite = fOriginalMaterials.cbegin();
            mimite != fOriginalMaterials.cend(); ++mimite )
  {
    fPhantomMaterials.push_back( (*mimite) );
  }
  //
  //---- Build and add new materials
  for( auto mppite= newMateDens.cbegin();
            mppite != newMateDens.cend(); ++mppite )
  {
    G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
    G4double saverdens = G4int(1000.001*averdens)/1000.;
    G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
           << averdens << " -> " << saverdens << " -> "
           << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000 << " "
           <<  G4int(1000*averdens)/1000. << G4endl;
 
    G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
           << averdens << " -> "  << saverdens << " -> "
           << G4UIcommand::ConvertToString(saverdens) << " Nvoxels "
    <<(*mppite).second->fNvoxels << G4endl;
    G4String mateName = ((*mppite).first).first->GetName() + "_" +
    G4UIcommand::ConvertToString(saverdens);
    fPhantomMaterials.push_back( BuildMaterialWithChangingDensity(
                                  (*mppite).first.first,
                                 G4float(averdens), mateName ) );
  }
}