DICOM/src/DicomPhantomZSliceHeader.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 //
0027 /// \file DicomPhantomZSliceHeader.cc
0028 /// \brief Implementation of the DicomPhantomZSliceHeader class
0029 //
0030
0031 #include "DicomPhantomZSliceHeader.hh"
0032
0033 #include "G4GeometryTolerance.hh"
0034 #include "G4LogicalVolume.hh"
0035 #include "G4Material.hh"
0036 #include "G4MaterialTable.hh"
0037 #include "G4NistManager.hh"
0038 #include "globals.hh"
0039
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
0041 DicomPhantomZSliceHeader::DicomPhantomZSliceHeader(const G4String& fname)
0042   : fNoVoxelsX(0),
0043     fNoVoxelsY(0),
0044     fNoVoxelsZ(0),
0045     fMinX(0),
0046     fMinY(0),
0047     fMinZ(0),
0048     fMaxX(0),
0049     fMaxY(0),
0050     fMaxZ(0),
0051     fFilename(fname),
0052     fSliceLocation(0)
0053 {}
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
0056 DicomPhantomZSliceHeader::~DicomPhantomZSliceHeader() {}
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
0059 DicomPhantomZSliceHeader::DicomPhantomZSliceHeader(std::ifstream& fin)
0060 {
0061   //----- Read material indices and names
0062   G4int nmate;
0063   G4String mateindex;
0064   G4String matename;
0065   fin >> nmate;
0066 #ifdef G4VERBOSE
0067   G4cout << " DicomPhantomZSliceHeader reading number of materials " << nmate << G4endl;
0068 #endif
0069
0070   for (G4int im = 0; im < nmate; im++) {
0071     fin >> mateindex >> matename;
0072 #ifdef G4VERBOSE
0073     // G4cout << " DicomPhantomZSliceHeader reading material "
0074     //  << im << " : "<< mateindex << "  " << matename << G4endl;
0075 #endif
0076
0077     if (!CheckMaterialExists(matename)) {
0078       G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader",
0079                   "A material is found in file that is not built in the C++ code",
0080                   FatalErrorInArgument, matename.c_str());
0081     }
0082
0083     fMaterialNames.push_back(matename);
0084   }
0085
0086   //----- Read number of voxels
0087   fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
0088 #ifdef G4VERBOSE
0089   G4cout << " Number of voxels " << fNoVoxelsX << " " << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
0090 #endif
0091
0092   //----- Read minimal and maximal extensions (= walls of phantom)
0093   fin >> fMinX >> fMaxX;
0094   fin >> fMinY >> fMaxY;
0095   fin >> fMinZ >> fMaxZ;
0096 #ifdef G4VERBOSE
0097   G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl << " Extension in Y " << fMinY
0098          << " " << fMaxY << G4endl << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
0099 #endif
0100
0101   fSliceLocation = 0.5 * (fMinZ + fMaxZ);
0102 }
0103
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0105 G4bool DicomPhantomZSliceHeader::CheckMaterialExists(const G4String& mateName)
0106 {
0107   const G4MaterialTable* matTab = G4Material::GetMaterialTable();
0108   std::vector<G4Material*>::const_iterator matite;
0109   for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
0110     if ((*matite)->GetName() == mateName) {
0111       return true;
0112     }
0113   }
0114
0115   G4Material* g4mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
0116   if (g4mate) {
0117     return false;
0118   }
0119   else {
0120     return true;
0121   }
0122 }
0123
0124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0125 void DicomPhantomZSliceHeader::operator+=(const DicomPhantomZSliceHeader& rhs)
0126 {
0127   *this = *this + rhs;
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
0131 DicomPhantomZSliceHeader DicomPhantomZSliceHeader::operator+(const DicomPhantomZSliceHeader& rhs)
0132 {
0133   //----- Check that both slices has the same dimensions
0134   if (fNoVoxelsX != rhs.GetNoVoxelsX() || fNoVoxelsY != rhs.GetNoVoxelsY()) {
0135     G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
0136         !!! Different number of voxels: "
0137            << "  X= " << fNoVoxelsX << " =? " << rhs.GetNoVoxelsX() << "  Y=  " << fNoVoxelsY
0138            << " =? " << rhs.GetNoVoxelsY() << "  Z=  " << fNoVoxelsZ << " =? " << rhs.GetNoVoxelsZ()
0139            << G4endl;
0140     G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader", "", FatalErrorInArgument, "");
0141   }
0142   //----- Check that both slices has the same extensions
0143   if (fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX() || fMinY != rhs.GetMinY()
0144       || fMaxY != rhs.GetMaxY())
0145   {
0146     G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
0147         !!! Different extensions: "
0148            << "  Xmin= " << fMinX << " =? " << rhs.GetMinX() << "  Xmax= " << fMaxX << " =? "
0149            << rhs.GetMaxX() << "  Ymin= " << fMinY << " =? " << rhs.GetMinY() << "  Ymax= " << fMaxY
0150            << " =? " << rhs.GetMaxY() << G4endl;
0151     G4Exception("DicomPhantomZSliceHeader::operator+", "", FatalErrorInArgument, "");
0152   }
0153
0154   //----- Check that both slices have the same materials
0155   std::vector<G4String> fMaterialNames2 = rhs.GetMaterialNames();
0156   if (fMaterialNames.size() != fMaterialNames2.size()) {
0157     G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
0158         !!! Different number of materials: "
0159            << fMaterialNames.size() << " =? " << fMaterialNames2.size() << G4endl;
0160     G4Exception("DicomPhantomZSliceHeader::operator+", "", FatalErrorInArgument, "");
0161   }
0162   for (unsigned int ii = 0; ii < fMaterialNames.size(); ii++) {
0163     if (fMaterialNames[ii] != fMaterialNames2[ii]) {
0164       G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:\
0165             !!! Different material number "
0166              << ii << " : " << fMaterialNames[ii] << " =? " << fMaterialNames2[ii] << G4endl;
0167       G4Exception("DicomPhantomZSliceHeader::operator+", "", FatalErrorInArgument, "");
0168     }
0169   }
0170
0171   //----- Check that the slices are contiguous in Z
0172   if (std::fabs(fMinZ - rhs.GetMaxZ()) > G4GeometryTolerance::GetInstance()->GetRadialTolerance()
0173       && std::fabs(fMaxZ - rhs.GetMinZ())
0174            > G4GeometryTolerance::GetInstance()->GetRadialTolerance())
0175   {
0176     G4cerr << "DicomPhantomZSliceHeader error adding two slice headers:!!!\
0177         Slices are not contiguous in Z "
0178            << "  Zmin= " << fMinZ << " & " << rhs.GetMinZ() << "  Zmax= " << fMaxZ << " & "
0179            << rhs.GetMaxZ() << G4endl;
0180     G4Exception("DicomPhantomZSliceHeader::operator+", "", FatalErrorInArgument, "");
0181   }
0182
0183   //----- Build slice header copying first one
0184   DicomPhantomZSliceHeader temp(*this);
0185
0186   //----- Add data from second slice header
0187   temp.SetMinZ(std::min(fMinZ, rhs.GetMinZ()));
0188   temp.SetMaxZ(std::max(fMaxZ, rhs.GetMaxZ()));
0189   temp.SetNoVoxelsZ(fNoVoxelsZ + rhs.GetNoVoxelsZ());
0190
0191   return temp;
0192 }
0193
0194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
0195 void DicomPhantomZSliceHeader::DumpToFile()
0196 {
0197   G4cout << "DicomPhantomZSliceHeader::Dumping Z Slice data to " << fFilename << "..." << G4endl;
0198   // sleep(5);
0199
0200   //  May seen counter-intuitive (dumping to file you are reading from), but
0201   //  the reason for this is modification slice spacing
0202   if (fMateIDs.size() == 0 || fValues.size() == 0) {
0203     ReadDataFromFile();
0204   }
0205
0206   std::ofstream out;
0207   out.open(fFilename.c_str());
0208
0209   if (!out) {
0210     G4String descript = "DicomPhantomZSliceHeader::DumpToFile: could not open " + fFilename;
0211     G4Exception(descript.c_str(), "", FatalException, "");
0212   }
0213
0214   out << fMaterialNames.size() << std::endl;
0215   for (unsigned int i = 0; i < fMaterialNames.size(); ++i) {
0216     out << i << " " << fMaterialNames.at(i) << std::endl;
0217   }
0218
0219   out << fNoVoxelsX << " " << fNoVoxelsY << " " << fNoVoxelsZ << std::endl;
0220   out << fMinX << " " << fMaxX << std::endl;
0221   out << fMinY << " " << fMaxY << std::endl;
0222   out << fMinZ << " " << fMaxZ << std::endl;
0223
0224   for (unsigned int i = 0; i < fMateIDs.size(); ++i) {
0225     Print(out, fMateIDs.at(i), " ");
0226   }
0227
0228   for (unsigned int i = 0; i < fValues.size(); ++i) {
0229     Print(out, fValues.at(i), " ", 6);
0230   }
0231
0232   out.close();
0233 }
0234
0235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0236 void DicomPhantomZSliceHeader::ReadDataFromFile()
0237 {
0238   std::ifstream in;
0239   in.open(fFilename.c_str());
0240
0241   if (!in) {
0242     G4String descript = "DicomPhantomZSliceHeader::DumpToFile: could not open " + fFilename;
0243     G4Exception(descript.c_str(), "", FatalException, "");
0244   }
0245
0246   G4int nMaterials;
0247   in >> nMaterials;
0248
0249   fMaterialNames.resize(nMaterials, "");
0250   for (G4int i = 0; i < nMaterials; ++i) {
0251     G4String str1, str2;
0252     in >> str1 >> str2;
0253     if (!IsInteger(str1)) {
0254       G4String descript = "String : " + str1 + " supposed to be integer";
0255       G4Exception(
0256         "DicomPhantomZSliceHeader::ReadDataFromFile - error in \
0257        formatting: missing material index",
0258         "", FatalException, descript.c_str());
0259     }
0260     G4int index = G4s2n<G4int>(str1);
0261     if (index > nMaterials || index < 0) {
0262       G4String descript = "Index : " + str1;
0263       G4Exception(
0264         "DicomPhantomZSliceHeader::ReadDataFromFile - error:\
0265             bad material index",
0266         "", FatalException, descript.c_str());
0267     }
0268     fMaterialNames[index] = str2;
0269   }
0270
0271   in >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
0272
0273   G4double tmpMinX, tmpMinY, tmpMinZ;
0274   G4double tmpMaxX, tmpMaxY, tmpMaxZ;
0275
0276   in >> tmpMinX >> tmpMaxX;
0277   in >> tmpMinY >> tmpMaxY;
0278   in >> tmpMinZ >> tmpMaxZ;
0279
0280   fMinX =
0281     (CheckConsistency(tmpMinX, fMinX, "Min X value")) ? fMinX : ((fMinX == 0) ? tmpMinX : fMinX);
0282   fMaxX =
0283     (CheckConsistency(tmpMaxX, fMaxX, "Max X value")) ? fMaxX : ((fMaxX == 0) ? tmpMaxX : fMaxX);
0284
0285   fMinY =
0286     (CheckConsistency(tmpMinY, fMinY, "Min Y value")) ? fMinY : ((fMinY == 0) ? tmpMinY : fMinY);
0287   fMaxY =
0288     (CheckConsistency(tmpMaxY, fMaxY, "Max Y value")) ? fMaxY : ((fMaxY == 0) ? tmpMaxY : fMaxY);
0289
0290   fMinZ =
0291     (CheckConsistency(tmpMinZ, fMinZ, "Min Z value")) ? fMinZ : ((fMinZ == 0) ? tmpMinZ : fMinZ);
0292   fMaxZ =
0293     (CheckConsistency(tmpMaxZ, fMaxZ, "Max Z value")) ? fMaxZ : ((fMaxZ == 0) ? tmpMaxZ : fMaxZ);
0294
0295   fMateIDs.clear();
0296   fValues.clear();
0297   fMateIDs.resize(fNoVoxelsY * fNoVoxelsZ, std::vector<G4int>(fNoVoxelsX, 0));
0298   fValues.resize(fNoVoxelsY * fNoVoxelsZ, std::vector<G4double>(fNoVoxelsX, 0.));
0299   for (G4int k = 0; k < fNoVoxelsZ; ++k) {
0300     for (G4int j = 0; j < fNoVoxelsY; ++j) {
0301       for (G4int i = 0; i < fNoVoxelsX; ++i) {
0302         G4int tmpMateID;
0303         in >> tmpMateID;
0304         G4int row = j * (k + 1);
0305         fMateIDs[row][i] = tmpMateID;
0306       }
0307     }
0308   }
0309
0310   for (G4int k = 0; k < fNoVoxelsZ; ++k) {
0311     for (G4int j = 0; j < fNoVoxelsY; ++j) {
0312       for (G4int i = 0; i < fNoVoxelsX; ++i) {
0313         G4double tmpValue;
0314         in >> tmpValue;
0315         G4int row = j * (k + 1);
0316         fValues[row][i] = tmpValue;
0317       }
0318     }
0319   }
0320
0321   in.close();
0322 }
0323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....