G4DicomReader/src/DicomVFileImage.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 #include "DicomVFileImage.hh"
0027
0028 #include "DicomFileStructure.hh"
0029 #include "DicomROI.hh"
0030 #include "dcmtk/dcmdata/dcdeftag.h"
0031 #include "dcmtk/dcmdata/dcfilefo.h"
0032 #include "dcmtk/dcmdata/dcpixel.h"
0033 #include "dcmtk/dcmdata/dcpixseq.h"
0034 #include "dcmtk/dcmdata/dcpxitem.h"
0035 #include "dcmtk/dcmrt/drtimage.h"
0036
0037 #include "G4GeometryTolerance.hh"
0038
0039 #include <set>
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042 DicomVFileImage::DicomVFileImage()
0043 {
0044   theFileMgr = DicomFileMgr::GetInstance();
0045 }
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 DicomVFileImage::DicomVFileImage(DcmDataset* dset) : DicomVFile(dset)
0049 {
0050   theFileMgr = DicomFileMgr::GetInstance();
0051 }
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054 void DicomVFileImage::ReadData()
0055 {
0056   std::vector<double> dImagePositionPatient = Read1Data(theDataset, DCM_ImagePositionPatient, 3);
0057   fLocation = dImagePositionPatient[2];
0058   std::vector<double> dSliceThickness = Read1Data(theDataset, DCM_SliceThickness, 1);
0059   std::vector<double> dPixelSpacing = Read1Data(theDataset, DCM_PixelSpacing, 2);
0060
0061   std::vector<double> dRows = Read1Data(theDataset, DCM_Rows, 1);
0062   std::vector<double> dColumns = Read1Data(theDataset, DCM_Columns, 1);
0063   fNoVoxelsY = dRows[0];
0064   fNoVoxelsX = dColumns[0];
0065   fNoVoxelsZ = 1;
0066
0067   fMinX = dImagePositionPatient[0];  // center of upper corner of pixel?
0068   fMaxX = dImagePositionPatient[0] + dColumns[0] * dPixelSpacing[0];
0069
0070   fMinY = dImagePositionPatient[1];
0071   fMaxY = dImagePositionPatient[1] + dRows[0] * dPixelSpacing[1];
0072
0073   fMinZ = dImagePositionPatient[2] - dSliceThickness[0] / 2.;
0074   fMaxZ = dImagePositionPatient[2] + dSliceThickness[0] / 2.;
0075   fVoxelDimX = dPixelSpacing[0];
0076   fVoxelDimY = dPixelSpacing[1];
0077   fVoxelDimZ = dSliceThickness[0];
0078
0079   if (DicomFileMgr::verbose >= debugVerb)
0080     G4cout << " DicomVFileImage::ReadData:  fNoVoxels " << fNoVoxelsX << " " << fNoVoxelsY << " "
0081            << fNoVoxelsZ << G4endl;
0082   if (DicomFileMgr::verbose >= debugVerb)
0083     G4cout << " DicomVFileImage::ReadData:  fMin " << fMinX << " " << fMinY << " " << fMinZ
0084            << G4endl;
0085   if (DicomFileMgr::verbose >= debugVerb)
0086     G4cout << " DicomVFileImage::ReadData:  fMax " << fMaxX << " " << fMaxY << " " << fMaxZ
0087            << G4endl;
0088   if (DicomFileMgr::verbose >= debugVerb)
0089     G4cout << " DicomVFileImage::ReadData:  fVoxelDim " << fVoxelDimX << " " << fVoxelDimY << " "
0090            << fVoxelDimZ << G4endl;
0091
0092   std::vector<double> dImageOrientationPatient =
0093     Read1Data(theDataset, DCM_ImageOrientationPatient, 6);
0094   fOrientationRows = G4ThreeVector(dImageOrientationPatient[0], dImageOrientationPatient[1],
0095                                    dImageOrientationPatient[2]);
0096   fOrientationColumns = G4ThreeVector(dImageOrientationPatient[3], dImageOrientationPatient[4],
0097                                       dImageOrientationPatient[5]);
0098
0099   if (fOrientationRows != G4ThreeVector(1, 0, 0) || fOrientationColumns != G4ThreeVector(0, 1, 0)) {
0100     G4cerr << " OrientationRows " << fOrientationRows << " OrientationColumns "
0101            << fOrientationColumns << G4endl;
0102     G4Exception("DicomVFileImage::ReadData", "DFCT0002", JustWarning,
0103                 "OrientationRows must be (1,0,0) and OrientationColumns (0,1,0), please contact "
0104                 "GAMOS authors");
0105   }
0106   fBitAllocated = Read1Data(theDataset, DCM_BitsAllocated, 1)[0];
0107   if (DicomFileMgr::verbose >= 4) G4cout << " BIT ALLOCATED " << fBitAllocated << G4endl;
0108
0109   std::vector<double> dRescaleSlope = Read1Data(theDataset, DCM_RescaleSlope, 1);
0110   if (dRescaleSlope.size() == 1) {
0111     fRescaleSlope = dRescaleSlope[0];
0112   }
0113   else {
0114     fRescaleSlope = 1;
0115   }
0116   std::vector<double> dRescaleIntercept = Read1Data(theDataset, DCM_RescaleIntercept, 1);
0117   if (dRescaleIntercept.size() == 1) {
0118     fRescaleIntercept = dRescaleIntercept[0];
0119   }
0120   else {
0121     fRescaleIntercept = 1;
0122   }
0123
0124   ReadPixelData();
0125 }
0126
0127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0128 void DicomVFileImage::ReadPixelData()
0129 {
0130   //  READING THE PIXELS :
0131   OFCondition result = EC_Normal;
0132   //---- CHECK IF DATA IS COMPRESSED
0133   DcmElement* element = NULL;
0134   result = theDataset->findAndGetElement(DCM_PixelData, element);
0135   if (result.bad() || element == NULL) {
0136     G4Exception("ReadData", "findAndGetElement(DCM_PixelData, ", FatalException,
0137                 ("Element PixelData not found: " + G4String(result.text())).c_str());
0138   }
0139   DcmPixelData* dpix = NULL;
0140   dpix = OFstatic_cast(DcmPixelData*, element);
0141   // If we have compressed data, we must utilize DcmPixelSequence
0142   //   in order to access it in raw format, e. g. for decompressing it
0143   //   with an external library.
0144   DcmPixelSequence* dseq = NULL;
0145   E_TransferSyntax xferSyntax = EXS_Unknown;
0146   const DcmRepresentationParameter* rep = NULL;
0147   // Find the key that is needed to access the right representation of the data within DCMTK
0148   dpix->getOriginalRepresentationKey(xferSyntax, rep);
0149   // Access original data representation and get result within pixel sequence
0150   result = dpix->getEncapsulatedRepresentation(xferSyntax, rep, dseq);
0151   if (result == EC_Normal)  // COMPRESSED DATA
0152   {
0153     G4Exception("DicomVFileImage::ReadData()", "DFCT004", FatalException,
0154                 "Compressed pixel data is not supported");
0155
0156     if (DicomFileMgr::verbose >= debugVerb)
0157       G4cout << " DicomVFileImage::ReadData:  result == EC_Normal Reading compressed data "
0158              << std::endl;
0159     DcmPixelItem* pixitem = NULL;
0160     // Access first frame (skipping offset table)
0161     for (int ii = 1; ii < 2; ii++) {
0162       OFCondition cond = dseq->getItem(pixitem, ii);
0163       if (!cond.good()) break;
0164       G4cout << ii << " PIX LENGTH " << pixitem->getLength() << G4endl;
0165     }
0166     if (pixitem == NULL) {
0167       G4Exception("ReadData", "dseq->getItem()", FatalException,
0168                   "No DcmPixelItem in DcmPixelSequence");
0169     }
0170     Uint8* pixData = NULL;
0171     // Get the length of this pixel item
0172     // (i.e. fragment, i.e. most of the time, the lenght of the frame)
0173     Uint32 length = pixitem->getLength();
0174     if (length == 0) {
0175       G4Exception("ReadData", "pixitem->getLength()", FatalException, "PixelData empty");
0176     }
0177
0178     if (DicomFileMgr::verbose >= debugVerb)
0179       G4cout << " DicomVFileImage::ReadData:  number of pixels " << length << G4endl;
0180     // Finally, get the compressed data for this pixel item
0181     result = pixitem->getUint8Array(pixData);
0182   }
0183   else {  // UNCOMPRESSED DATA
0184     if (fBitAllocated == 8) {  // Case 8 bits :
0185       Uint8* pixData = NULL;
0186       if (!(element->getUint8Array(pixData)).good()) {
0187         G4Exception("ReadData", "getUint8Array pixData, ", FatalException,
0188                     ("PixelData not found: " + G4String(result.text())).c_str());
0189       }
0190       for (int ir = 0; ir < fNoVoxelsY; ir++) {
0191         for (int ic = 0; ic < fNoVoxelsX; ic++) {
0192           fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
0193         }
0194       }
0195     }
0196     else if (fBitAllocated == 16) {  // Case 16 bits :
0197       Uint16* pixData = NULL;
0198       if (!(element->getUint16Array(pixData)).good()) {
0199         G4Exception("ReadData", "getUint16Array pixData, ", FatalException,
0200                     ("PixelData not found: " + G4String(result.text())).c_str());
0201       }
0202       for (int ir = 0; ir < fNoVoxelsY; ir++) {
0203         for (int ic = 0; ic < fNoVoxelsX; ic++) {
0204           fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
0205         }
0206       }
0207     }
0208     else if (fBitAllocated == 32) {  // Case 32 bits :
0209       Uint32* pixData = NULL;
0210       if (!(element->getUint32Array(pixData)).good()) {
0211         G4Exception("ReadData", "getUint32Array pixData, ", FatalException,
0212                     ("PixelData not found: " + G4String(result.text())).c_str());
0213       }
0214       for (int ir = 0; ir < fNoVoxelsY; ir++) {
0215         for (int ic = 0; ic < fNoVoxelsX; ic++) {
0216           fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
0217         }
0218       }
0219     }
0220   }
0221 }
0222
0223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0224 void DicomVFileImage::operator+=(const DicomVFileImage& rhs)
0225 {
0226   *this = *this + rhs;
0227 }
0228
0229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0230 DicomVFileImage DicomVFileImage::operator+(const DicomVFileImage& rhs)
0231 {
0232   //----- Check that both slices has the same dimensions
0233   if (fNoVoxelsX != rhs.GetNoVoxelsX() || fNoVoxelsY != rhs.GetNoVoxelsY()) {
0234     G4cerr << "DicomVFileImage error adding two slice headers:\
0235         !!! Different number of voxels: "
0236            << "  X= " << fNoVoxelsX << " =? " << rhs.GetNoVoxelsX() << "  Y=  " << fNoVoxelsY
0237            << " =? " << rhs.GetNoVoxelsY() << "  Z=  " << fNoVoxelsZ << " =? " << rhs.GetNoVoxelsZ()
0238            << G4endl;
0239     G4Exception("DicomVFileImage::DicomVFileImage", "", FatalErrorInArgument, "");
0240   }
0241   //----- Check that both slices has the same extensions
0242   if (fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX() || fMinY != rhs.GetMinY()
0243       || fMaxY != rhs.GetMaxY())
0244   {
0245     G4cerr << "DicomVFileImage error adding two slice headers:\
0246         !!! Different extensions: "
0247            << "  Xmin= " << fMinX << " =? " << rhs.GetMinX() << "  Xmax= " << fMaxX << " =? "
0248            << rhs.GetMaxX() << "  Ymin= " << fMinY << " =? " << rhs.GetMinY() << "  Ymax= " << fMaxY
0249            << " =? " << rhs.GetMaxY() << G4endl;
0250     G4Exception("DicomVFileImage::operator+", "", FatalErrorInArgument, "");
0251   }
0252
0253   //----- Check that both slices has the same orientations
0254   if (fOrientationRows != rhs.GetOrientationRows()
0255       || fOrientationColumns != rhs.GetOrientationColumns())
0256   {
0257     G4cerr << "DicomVFileImage error adding two slice headers: !!!\
0258         Slices have different orientations "
0259            << "  Orientation Rows = " << fOrientationRows << " & " << rhs.GetOrientationRows()
0260            << "  Orientation Columns " << fOrientationColumns << " & "
0261            << rhs.GetOrientationColumns() << G4endl;
0262     G4Exception("DicomVFileImage::operator+", "", FatalErrorInArgument, "");
0263   }
0264
0265   //----- Check that the slices are contiguous in Z
0266   if (std::fabs(fMinZ - rhs.GetMaxZ()) > G4GeometryTolerance::GetInstance()->GetRadialTolerance()
0267       && std::fabs(fMaxZ - rhs.GetMinZ())
0268            > G4GeometryTolerance::GetInstance()->GetRadialTolerance())
0269   {
0270     G4cerr << "DicomVFileImage error adding two slice headers: !!!\
0271         Slices are not contiguous in Z "
0272            << "  Zmin= " << fMinZ << " & " << rhs.GetMinZ() << "  Zmax= " << fMaxZ << " & "
0273            << rhs.GetMaxZ() << G4endl;
0274     G4Exception("DicomVFileImage::operator+", "", JustWarning, "");
0275   }
0276
0277   //----- Build slice header copying first one
0278   DicomVFileImage temp(*this);
0279
0280   //----- Add data from second slice header
0281   temp.SetMinZ(std::min(fMinZ, rhs.GetMinZ()));
0282   temp.SetMaxZ(std::max(fMaxZ, rhs.GetMaxZ()));
0283   temp.SetNoVoxelsZ(fNoVoxelsZ + rhs.GetNoVoxelsZ());
0284
0285   return temp;
0286 }
0287
0288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0289 void DicomVFileImage::DumpHeaderToTextFile(std::ofstream& fout)
0290 {
0291   if (DicomFileMgr::verbose >= warningVerb)
0292     G4cout << fLocation << " DumpHeaderToTextFile " << fFileName << " " << fHounsfieldV.size()
0293            << G4endl;
0294
0295   G4String fName = fFileName.substr(0, fFileName.length() - 3) + "g4dcm";
0296   std::ofstream out(fName.c_str());
0297
0298   if (DicomFileMgr::verbose >= warningVerb)
0299     G4cout << "### DicomVFileImage::Dumping Z Slice header to Text file " << G4endl;
0300
0301   G4int fCompress = theFileMgr->GetCompression();
0302   fout << fNoVoxelsX / fCompress << " " << fNoVoxelsY / fCompress << " " << fNoVoxelsZ << std::endl;
0303   fout << fMinX << " " << fMaxX << std::endl;
0304   fout << fMinY << " " << fMaxY << std::endl;
0305   fout << fMinZ << " " << fMaxZ << std::endl;
0306 }
0307
0308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0309 void DicomVFileImage::Print(std::ostream& out)
0310 {
0311   G4int fCompress = theFileMgr->GetCompression();
0312   out << "@@ CT Slice " << fLocation << G4endl;
0313
0314   out << "@ NoVoxels " << fNoVoxelsX / fCompress << " " << fNoVoxelsY / fCompress << " "
0315       << fNoVoxelsZ << G4endl;
0316   out << "@ DIM X: " << fMinX << " " << fMaxX << " Y: " << fMinY << " " << fMaxY << " Z: " << fMinZ
0317       << " " << fMaxZ << G4endl;
0318 }