DICOM/src/DicomNestedPhantomParameterisation.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file medical/DICOM/src/DicomNestedPhantomParameterisation.cc
0027 /// \brief Implementation of the DicomNestedPhantomParameterisation class
0028 //
0029 //
0030
0031 #include "DicomNestedPhantomParameterisation.hh"
0032
0033 #include "DicomHandler.hh"
0034
0035 #include "G4Box.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4Material.hh"
0038 #include "G4ThreeVector.hh"
0039 #include "G4VPhysicalVolume.hh"
0040 #include "G4VTouchable.hh"
0041 #include "G4VVisManager.hh"
0042 #include "G4VisAttributes.hh"
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045
0046 G4String DicomNestedPhantomParameterisation::fDefaultColorFile =
0047   DicomHandler::GetDicomDataPath() + "/ColourMap.dat";
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050 DicomNestedPhantomParameterisation::DicomNestedPhantomParameterisation(
0051   const G4ThreeVector& voxelSize, std::vector<G4Material*>& mat, G4int fnZ_, G4int fnY_, G4int fnX_,
0052   G4String colorfile)
0053   :  // G4VNestedParameterisation(),
0054     fdX(voxelSize.x()),
0055     fdY(voxelSize.y()),
0056     fdZ(voxelSize.z()),
0057     fnX(fnX_),
0058     fnY(fnY_),
0059     fnZ(fnZ_),
0060     fMaterials(mat),
0061     fMaterialIndices(0)
0062 {
0063   ReadColourData(colorfile);
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067 DicomNestedPhantomParameterisation::~DicomNestedPhantomParameterisation() {}
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070 void DicomNestedPhantomParameterisation::ReadColourData(G4String colourFile)
0071 {
0072   //----- Add a G4VisAttributes for materials not defined in file;
0073   G4VisAttributes* blankAtt = new G4VisAttributes;
0074   blankAtt->SetVisibility(FALSE);
0075   fColours["Default"] = blankAtt;
0076
0077   std::ifstream fin(colourFile.c_str());
0078   G4int nMate;
0079   G4String mateName;
0080   G4double cred, cgreen, cblue, copacity;
0081   fin >> nMate;
0082   for (G4int ii = 0; ii < nMate; ii++) {
0083     fin >> mateName;
0084     if (fin.eof()) break;
0085     fin >> cred >> cgreen >> cblue >> copacity;
0086     G4Colour colour(cred, cgreen, cblue, copacity);
0087     G4VisAttributes* visAtt = new G4VisAttributes(colour);
0088     visAtt->SetVisibility(true);
0089     fColours[mateName] = visAtt;
0090     fColours2[ii] = new G4VisAttributes(*visAtt);
0091   }
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095 void DicomNestedPhantomParameterisation::SetNoVoxels(unsigned int nx, unsigned int ny,
0096                                                      unsigned int nz)
0097 {
0098   fnX = nx;
0099   fnY = ny;
0100   fnZ = nz;
0101 }
0102
0103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0104 G4Material* DicomNestedPhantomParameterisation::ComputeMaterial(G4VPhysicalVolume* physVol,
0105                                                                 const G4int iz,
0106                                                                 const G4VTouchable* parentTouch)
0107 {
0108   // protection for initialization and vis at idle state
0109   //
0110   if (parentTouch == nullptr) return fMaterials[0];
0111
0112   // Copy number of voxels.
0113   // Copy number of X and Y are obtained from replication number.
0114   // Copy nymber of Z is the copy number of current voxel.
0115   G4int ix = parentTouch->GetReplicaNumber(0);
0116   G4int iy = parentTouch->GetReplicaNumber(1);
0117
0118   G4int copyID = ix + fnX * iy + fnX * fnY * iz;
0119
0120   std::size_t matIndex = GetMaterialIndex(copyID);
0121   G4Material* mate = fMaterials[matIndex];
0122
0123   if (G4VVisManager::GetConcreteInstance() && physVol) {
0124     G4String mateName = fMaterials.at(matIndex)->GetName();
0125     std::string::size_type iuu = mateName.find("__");
0126     if (iuu != std::string::npos) mateName = mateName.substr(0, iuu);
0127
0128     if (0 < fColours.count(mateName))
0129       physVol->GetLogicalVolume()->SetVisAttributes(fColours.find(mateName)->second);
0130     else {
0131       bool found = false;
0132       for (const auto& itr : fColours) {
0133         G4String mat_color = itr.first;
0134         auto len = mat_color.length();
0135         if (mateName.find(mat_color) == 0 && mateName.length() > len && mateName[len] == '_') {
0136           physVol->GetLogicalVolume()->SetVisAttributes(fColours.find(mat_color)->second);
0137           found = true;
0138         }
0139         if (found) break;
0140       }
0141       if (!found) {
0142         static uintmax_t n = 0;
0143         if (n++ < 100)
0144           G4cout << "Unknown material name " << mateName << " for index " << matIndex << G4endl;
0145         if (fColours2.find(matIndex) != fColours2.end())
0146           physVol->GetLogicalVolume()->SetVisAttributes(fColours2.find(matIndex)->second);
0147         else
0148           physVol->GetLogicalVolume()->SetVisAttributes(fColours.begin()->second);
0149       }
0150     }
0151     physVol->GetLogicalVolume()->SetMaterial(mate);
0152   }
0153
0154   return mate;
0155 }
0156
0157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0158 unsigned int DicomNestedPhantomParameterisation::GetMaterialIndex(unsigned int copyNo) const
0159 {
0160   // return *(fMaterialIndices+copyNo);
0161   return unsigned(fMaterialIndices[copyNo]);
0162 }
0163
0164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0165 // Number of Materials
0166 // Material scanner is required for preparing physics tables and so on before
0167 // starting simulation, so that G4 has to know number of materials.
0168 //
0169 G4int DicomNestedPhantomParameterisation::GetNumberOfMaterials() const
0170 {
0171   return G4int(fMaterials.size());
0172 }
0173
0174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0175 //
0176 // GetMaterial
0177 //  This is needed for material scanner and realizing geometry.
0178 //
0179 G4Material* DicomNestedPhantomParameterisation::GetMaterial(G4int i) const
0180 {
0181   return fMaterials[i];
0182 }
0183
0184 //
0185 // Transformation of voxels.
0186 //
0187 void DicomNestedPhantomParameterisation::ComputeTransformation(const G4int copyNo,
0188                                                                G4VPhysicalVolume* physVol) const
0189 {
0190   // Position of voxels.
0191   // x and y positions are already defined in DetectorConstruction by using
0192   // replicated volume. Here only we need to define is z positions of voxels.
0193   physVol->SetTranslation(
0194     G4ThreeVector(0., 0., (2. * static_cast<double>(copyNo) + 1.) * fdZ - fdZ * fnZ));
0195 }
0196
0197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0198 //
0199 // Dimensions are always same in this RE02 example.
0200 //
0201 void DicomNestedPhantomParameterisation::ComputeDimensions(G4Box& box, const G4int,
0202                                                            const G4VPhysicalVolume*) const
0203 {
0204   box.SetXHalfLength(fdX);
0205   box.SetYHalfLength(fdY);
0206   box.SetZHalfLength(fdZ);
0207 }