DicomHandler Class Reference

#include <Doxymodules_medical.h>

Public Member Functions
	~DicomHandler ()
G4int	ReadFile (FILE *, char *)
G4int	ReadData (FILE *, char *)
void	CheckFileFormat ()

Static Public Member Functions
static DicomHandler *	Instance ()
static G4String	GetDicomDataPath ()
static G4String	GetDicomDataFile ()

Private Member Functions
	DicomHandler ()
template<class Type >
void	GetValue (char *, Type &)
void	ReadCalibration ()
void	GetInformation (G4int &, char *)
G4float	Pixel2density (G4int pixel)
void	ReadMaterialIndices (std::ifstream &finData)
unsigned int	GetMaterialIndex (G4float density)
void	StoreData (std::ofstream &foutG4DCM)
void	StoreData (DicomPhantomZSliceHeader *dcmPZSH)
G4int	read_defined_nested (FILE *, G4int)
void	read_undefined_nested (FILE *)
void	read_undefined_item (FILE *)

Private Attributes
const G4int	DATABUFFSIZE
const G4int	LINEBUFFSIZE
const G4int	FILENAMESIZE
short	fCompression
G4int	fNFiles
short	fRows
short	fColumns
short	fBitAllocated
G4int	fMaxPixelValue
G4int	fMinPixelValue
G4double	fPixelSpacingX
G4double	fPixelSpacingY
G4double	fSliceThickness
G4double	fSliceLocation
G4int	fRescaleIntercept
G4int	fRescaleSlope
G4bool	fLittleEndian
G4bool	fImplicitEndian
short	fPixelRepresentation
G4int **	fTab
std::map< G4float, G4String >	fMaterialIndices
G4int	fNbrequali
G4double *	fValueDensity
G4double *	fValueCT
G4bool	fReadCalibration
DicomPhantomZSliceMerged *	fMergedSlices
G4String	fDriverFile
G4String	fCt2DensityFile

Static Private Attributes
static DicomHandler *	fInstance = 0
	DicomHandler.cc :

Detailed Description

Definition at line 21 of file Doxymodules_medical.h.

Constructor & Destructor Documentation

~DicomHandler()

DicomHandler::~DicomHandler

(

)

Definition at line 164 of file DicomHandler.cc.

{}

DicomHandler()

DicomHandler::DicomHandler ( )

private

Definition at line 148 of file DicomHandler.cc.

:DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
 fCompression(0), fNFiles(0), fRows(0), fColumns(0),
 fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
 fPixelSpacingX(0.), fPixelSpacingY(0.),fSliceThickness(0.), 
 fSliceLocation(0.), fRescaleIntercept(0), fRescaleSlope(0),
 fLittleEndian(true),fImplicitEndian(false),fPixelRepresentation(0),
 fNbrequali(0),fValueDensity(NULL),fValueCT(NULL),
 fReadCalibration(false),fMergedSlices(NULL),
 fDriverFile("Data.dat"),fCt2DensityFile("CT2Density.dat")
{
 fMergedSlices = new DicomPhantomZSliceMerged;
}

Member Function Documentation

Instance()

DicomHandler * DicomHandler::Instance ( )

static

Definition at line 71 of file DicomHandler.cc.

{
  if (fInstance == 0)
  {
    static DicomHandler dicomhandler;
    fInstance = &dicomhandler;
  }
  return fInstance;
}

ReadFile()

G4int DicomHandler::ReadFile	(	FILE *	dicom,
		char *	filename2
	)

Definition at line 169 of file DicomHandler.cc.

{
    G4cout << " ReadFile " << filename2 << G4endl;
 
    G4int returnvalue = 0; size_t rflag = 0;
    char * buffer = new char[LINEBUFFSIZE];
 
    fImplicitEndian = false;
    fLittleEndian = true;
 
    rflag = std::fread( buffer, 1, 128, dicom ); // The first 128 bytes
                                                 //are not important
                                                 // Reads the "DICOM" letters
    rflag = std::fread( buffer, 1, 4, dicom );
    // if there is no preamble, the FILE pointer is rewinded.
    if(std::strncmp("DICM", buffer, 4) != 0) {
        std::fseek(dicom, 0, SEEK_SET);
        fImplicitEndian = true;
    }
 
    short readGroupId;    // identify the kind of input data
    short readElementId;  // identify a particular type information
    short elementLength2; // deal with element length in 2 bytes
                          //unsigned int elementLength4; 
    // deal with element length in 4 bytes
    unsigned long elementLength4; // deal with element length in 4 bytes
 
    char * data = new char[DATABUFFSIZE];
 
    // Read information up to the pixel data
    while(true) {
      
      //Reading groups and elements :
      readGroupId = 0;
      readElementId = 0;
      // group ID
      rflag = std::fread(buffer, 2, 1, dicom);
      GetValue(buffer, readGroupId);
      // element ID
      rflag = std::fread(buffer, 2, 1, dicom);
      GetValue(buffer, readElementId);
      
      // Creating a tag to be identified afterward
      G4int tagDictionary = readGroupId*0x10000 + readElementId;
      
      // beginning of the pixels
      if(tagDictionary == 0x7FE00010) {
        // Following 2 fread's are modifications to original DICOM example
        // (Jonathan Madsen)
        if(!fImplicitEndian)
            rflag = std::fread(buffer,2,1,dicom);   // Reserved 2 bytes
        // (not used for pixels)
        rflag = std::fread(buffer,4,1,dicom);   // Element Length  
        // (not used for pixels)
        break;      // Exit to ReadImageData()
      }
      
      // VR or element length
      rflag = std::fread(buffer,2,1,dicom);
      GetValue(buffer, elementLength2);
      
      // If value representation (VR) is OB, OW, SQ, UN, added OF and UT
      //the next length is 32 bits
      if((elementLength2 == 0x424f ||     // "OB"
        elementLength2 == 0x574f ||     // "OW"
        elementLength2 == 0x464f ||     // "OF"
        elementLength2 == 0x5455 ||     // "UT"
        elementLength2 == 0x5153 ||     // "SQ"
        elementLength2 == 0x4e55) &&    // "UN"
       !fImplicitEndian ) {             // explicit VR
      
      rflag = std::fread(buffer, 2, 1, dicom); // Skip 2 reserved bytes
      
      // element length
      rflag = std::fread(buffer, 4, 1, dicom);
      GetValue(buffer, elementLength4);
      
      if(elementLength2 == 0x5153)
        {
          if(elementLength4 == 0xFFFFFFFF)
            {
            read_undefined_nested( dicom );
            elementLength4=0;
            }  else{
            if(read_defined_nested( dicom, elementLength4 )==0){
            G4Exception("DicomHandler::ReadFile",
                      "DICOM001",
                      FatalException,
                      "Function read_defined_nested() failed!");
            }
          }
        } else  {
        // Reading the information with data
        rflag = std::fread(data, elementLength4,1,dicom);
      }
      
      }  else {
 
      //  explicit with VR different than previous ones
      if(!fImplicitEndian || readGroupId == 2) {
        
        //G4cout << "Reading  DICOM files with Explicit VR"<< G4endl;
        // element length (2 bytes)
        rflag = std::fread(buffer, 2, 1, dicom);
        GetValue(buffer, elementLength2);
        elementLength4 = elementLength2;
        
        rflag = std::fread(data, elementLength4, 1, dicom);
        
      } else {                                  // Implicit VR
        
        //G4cout << "Reading  DICOM files with Implicit VR"<< G4endl;
        
        // element length (4 bytes)
        if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
          G4Exception("DicomHandler::ReadFile",
                  "DICOM001",
                  FatalException,
                  "fseek failed");
        }
    
        rflag = std::fread(buffer, 4, 1, dicom);
        GetValue(buffer, elementLength4);
        
        //G4cout <<  std::hex<< elementLength4 << G4endl;
        
        if(elementLength4 == 0xFFFFFFFF)
          {
            read_undefined_nested(dicom);
            elementLength4=0;
          }  else{
          rflag = std::fread(data, elementLength4, 1, dicom);
        }
        
      }
      }
      
      // NULL termination
      data[elementLength4] = '\0';
      
      // analyzing information
      GetInformation(tagDictionary, data);
    }
    
    G4String fnameG4DCM = G4String(filename2) + ".g4dcm";
 
    // Perform functions originally written straight to file
    DicomPhantomZSliceHeader* zslice = new DicomPhantomZSliceHeader(fnameG4DCM);
    for(auto ite = fMaterialIndices.cbegin();
             ite != fMaterialIndices.cend(); ++ite)
    {
      zslice->AddMaterial(ite->second);
    }
    
    zslice->SetNoVoxelsX(fColumns/fCompression);
    zslice->SetNoVoxelsY(fRows/fCompression);
    zslice->SetNoVoxelsZ(1);
 
    zslice->SetMinX(-fPixelSpacingX*fColumns/2.);
    zslice->SetMaxX(fPixelSpacingX*fColumns/2.);
 
    zslice->SetMinY(-fPixelSpacingY*fRows/2.);
    zslice->SetMaxY(fPixelSpacingY*fRows/2.);
 
    zslice->SetMinZ(fSliceLocation-fSliceThickness/2.);
    zslice->SetMaxZ(fSliceLocation+fSliceThickness/2.);
 
    //===
 
    ReadData( dicom, filename2 );
 
    // DEPRECIATED
    //StoreData( foutG4DCM );
    //foutG4DCM.close();
 
    StoreData( zslice );
 
    // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency
    //zslice->DumpToFile();
 
    fMergedSlices->AddZSlice(zslice);
 
    //
    delete [] buffer;
    delete [] data;
 
    if (rflag) return returnvalue;
    return returnvalue;
}

ReadData()

G4int DicomHandler::ReadData	(	FILE *	dicom,
		char *	filename2
	)

Definition at line 679 of file DicomHandler.cc.

{
  G4int returnvalue = 0; size_t rflag = 0;
  
  //  READING THE PIXELS
  
  fTab = new G4int*[fRows];
  for ( G4int i = 0; i < fRows; ++i )
  {
    fTab[i] = new G4int[fColumns];
  }
  
  if(fBitAllocated == 8) { // Case 8 bits :
    
    std::printf("@@@ Error! Picture != 16 bits...\n");
    std::printf("@@@ Error! Picture != 16 bits...\n");
    std::printf("@@@ Error! Picture != 16 bits...\n");
    
    unsigned char ch = 0;
    
    for(G4int j = 0; j < fRows; ++j) {
      for(G4int i = 0; i < fColumns; ++i) {
        rflag = std::fread( &ch, 1, 1, dicom);
        fTab[j][i] = ch*fRescaleSlope + fRescaleIntercept;
      }
    }
    returnvalue = 1;
    
  } else { //  from 12 to 16 bits :
    char sbuff[2];
    short pixel;
    for( G4int j = 0; j < fRows; ++j) {
      for( G4int i = 0; i < fColumns; ++i) {
        rflag = std::fread(sbuff, 2, 1, dicom);
        GetValue(sbuff, pixel);
        fTab[j][i] = pixel*fRescaleSlope + fRescaleIntercept;
      }
    }
  }
  
  // Creation of .g4 files wich contains averaged density data
  G4String nameProcessed = filename2 + G4String(".g4dcmb");
  FILE* fileOut = std::fopen(nameProcessed.c_str(),"w+b");
  
  G4cout << "### Writing of "<< nameProcessed << " ###\n";
 
  unsigned int nMate = fMaterialIndices.size();
  rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
  //--- Write materials
  for( auto ite = fMaterialIndices.cbegin();
            ite != fMaterialIndices.cend(); ++ite )
  {
    G4String mateName = (*ite).second;
    for( std::size_t ii = (*ite).second.length(); ii < 40; ++ii )
    {
      mateName += " ";
    }         //mateName = const_cast<char*>(((*ite).second).c_str());
    
    const char* mateNameC = mateName.c_str();
    rflag = std::fwrite(mateNameC, sizeof(char),40, fileOut);
  }
 
  unsigned int fRowsC = fRows/fCompression;
  unsigned int fColumnsC = fColumns/fCompression;
  unsigned int planesC = 1;
  G4float pixelLocationXM = -G4float(fPixelSpacingX*fColumns/2.);
  G4float pixelLocationXP = G4float(fPixelSpacingX*fColumns/2.);
  G4float pixelLocationYM = -G4float(fPixelSpacingY*fRows/2.);
  G4float pixelLocationYP = G4float(fPixelSpacingY*fRows/2.);
  G4float fSliceLocationZM = G4float(fSliceLocation-fSliceThickness/2.);
  G4float fSliceLocationZP = G4float(fSliceLocation+fSliceThickness/2.);
  //--- Write number of voxels (assume only one voxel in Z)
  rflag = std::fwrite(&fRowsC, sizeof(unsigned int), 1, fileOut);
  rflag = std::fwrite(&fColumnsC, sizeof(unsigned int), 1, fileOut);
  rflag = std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut);
  //--- Write minimum and maximum extensions
  rflag = std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut);
  rflag = std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut);
  rflag = std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut);
  rflag = std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut);
  rflag = std::fwrite(&fSliceLocationZM, sizeof(G4float), 1, fileOut);
  rflag = std::fwrite(&fSliceLocationZP, sizeof(G4float), 1, fileOut);
  // rflag = std::fwrite(&fCompression, sizeof(unsigned int), 1, fileOut);
  
  std::printf("%8i   %8i\n",fRows,fColumns);
  std::printf("%8f   %8f\n",fPixelSpacingX,fPixelSpacingY);
  std::printf("%8f\n", fSliceThickness);
  std::printf("%8f\n", fSliceLocation);
  std::printf("%8i\n", fCompression);
  
  G4int compSize = fCompression;
  G4int mean;
  G4float density;
  G4bool overflow = false;
  
  //----- Write index of material for each pixel
  if(compSize == 1) { // no fCompression: each pixel has a density value)
    for( G4int ww = 0; ww < fRows; ++ww) {
      for( G4int xx = 0; xx < fColumns; ++xx) {
        mean = fTab[ww][xx];
        density = Pixel2density(mean);
        unsigned int mateID = GetMaterialIndex( density );
        rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
      }
    }
    
  } else {
    // density value is the average of a square region of
    // fCompression*fCompression pixels
    for(G4int ww = 0; ww < fRows ;ww += compSize ) {
      for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
        overflow = false;
        mean = 0;
        for(G4int sumx = 0; sumx < compSize; ++sumx) {
          for(G4int sumy = 0; sumy < compSize; ++sumy) {
            if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
            mean += fTab[ww+sumy][xx+sumx];
          }
          if(overflow) break;
        }
        mean /= compSize*compSize;
        
        if(!overflow) {
          density = Pixel2density(mean);
          unsigned int mateID = GetMaterialIndex( density );
          rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
        }
      }
      
    }
  }
  
  //----- Write density for each pixel
  if(compSize == 1) { // no fCompression: each pixel has a density value)
    for( G4int ww = 0; ww < fRows; ++ww) {
      for( G4int xx = 0; xx < fColumns; ++xx) {
        mean = fTab[ww][xx];
        density = Pixel2density(mean);
        rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
      }
    }
    
  } else {
    // density value is the average of a square region of
    // fCompression*fCompression pixels
    for(G4int ww = 0; ww < fRows ;ww += compSize ) {
      for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
        overflow = false;
        mean = 0;
        for(G4int sumx = 0; sumx < compSize; ++sumx) {
          for(G4int sumy = 0; sumy < compSize; ++sumy) {
            if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
            mean += fTab[ww+sumy][xx+sumx];
          }
          if(overflow) break;
        }
        mean /= compSize*compSize;
        
        if(!overflow) {
          density = Pixel2density(mean);
          rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
        }
      }
      
    }
  }
  
  rflag = std::fclose(fileOut);
  
  delete [] nameProcessed;
 
  /*    for ( G4int i = 0; i < fRows; i ++ ) {
        delete [] fTab[i];
        }
     delete [] fTab;
  */
  
  if (rflag) return returnvalue;
  return returnvalue;
}

CheckFileFormat()

void DicomHandler::CheckFileFormat

(

)

Definition at line 962 of file DicomHandler.cc.

{
  std::ifstream checkData(fDriverFile.c_str());
  char * oneLine = new char[128];
  
  if(!(checkData.is_open())) { //Check existance of Data.dat
    
    G4String message = 
      "\nDicomG4 needs Data.dat (or another driver file specified";
    message += " in command line):\n";
    message += "\tFirst line: number of image pixel for a voxel (G4Box)\n";
    message += "\tSecond line: number of images (CT slices) to read\n";
    message += "\tEach following line contains the name of a Dicom image";
    message += " except for the .dcm extension";
    G4Exception("DicomHandler::ReadFile",
                "DICOM001",
                FatalException,
                message.c_str());
  }
  
  checkData >> fCompression;
  checkData >> fNFiles;
  G4String oneName;
  checkData.getline(oneLine,100);
  std::ifstream testExistence;
  G4bool existAlready = true;
  for(G4int rep = 0; rep < fNFiles; ++rep) {
    checkData.getline(oneLine,100);
    oneName = oneLine;
    oneName += ".g4dcm"; // create dicomFile.g4dcm
    G4cout << fNFiles << " test file " << oneName << G4endl;
    testExistence.open(oneName.data());
    if(!(testExistence.is_open())) {
      existAlready = false;
      testExistence.clear();
      testExistence.close();
    }
    testExistence.clear();
    testExistence.close();
  }
  
  ReadMaterialIndices( checkData );
  
  checkData.close();
  delete [] oneLine;
  
  if( existAlready == false  ) { // The files *.g4dcm have to be created
    
    G4cout << "\nAll the necessary images were not found in processed form "
           << ", starting with .dcm images\n";
    
    FILE * dicom;
    FILE * lecturePref;
    char * fCompressionc = new char[LINEBUFFSIZE];
    char * maxc = new char[LINEBUFFSIZE];
    //char name[300], inputFile[300];
    char * inputFile = new char[FILENAMESIZE];
    G4int rflag;
    lecturePref = std::fopen(fDriverFile.c_str(),"r");
 
    rflag = std::fscanf(lecturePref,"%s",fCompressionc);
    fCompression = atoi(fCompressionc);
    rflag = std::fscanf(lecturePref,"%s",maxc);
    fNFiles = atoi(maxc);
    G4cout << " fNFiles " << fNFiles << G4endl;
  
///////////////////// 
 
#ifdef DICOM_USE_HEAD
    for( G4int i = 1; i <= fNFiles; ++i ) 
    { // Begin loop on filenames
     rflag = std::fscanf(lecturePref,"%s",inputFile);
      G4String path = GetDicomDataPath() + "/";
    
     G4String name = inputFile + G4String(".dcm");
     //Writes the results to a character string buffer.
     
     G4String name2 = path + name;
     //  Open input file and give it to gestion_dicom :
     G4cout << "### Opening " << name2 << " and reading :\n";
     dicom = std::fopen(name2.c_str(),"rb");
     // Reading the .dcm in two steps:
     //      1.  reading the header
     //      2. reading the pixel data and store the density in Moyenne.dat
     if( dicom != 0 ) {
        ReadFile(dicom,inputFile);
      } else {
        G4cout << "\nError opening file : " << name2 << G4endl;
      }
      rflag = std::fclose(dicom);
    }
#else
 
     for( G4int i = 1; i <= fNFiles; ++i ) 
     { // Begin loop on filenames
      rflag = std::fscanf(lecturePref,"%s",inputFile);
 
      G4String name = inputFile + G4String(".dcm");
      //Writes the results to a character string buffer.
      
       //G4cout << "check: " << name << G4endl;
      //  Open input file and give it to gestion_dicom :
      G4cout << "### Opening " << name << " and reading :\n";
      dicom = std::fopen(name.c_str(),"rb");
      // Reading the .dcm in two steps:
      //      1.  reading the header
      //      2. reading the pixel data and store the density in Moyenne.dat
      if( dicom != 0 ) {
        ReadFile(dicom,inputFile);
      } else {
        G4cout << "\nError opening file : " << name << G4endl;
      }
      rflag = std::fclose(dicom);
    }
#endif
 
    rflag = std::fclose(lecturePref);
    
    // Checks the spacing is correct. Dumps to files
    fMergedSlices->CheckSlices();
    
    delete [] fCompressionc;
    delete [] maxc;
    delete [] inputFile;
    if (rflag) return;
    
  }
  
  if(fValueDensity) { delete [] fValueDensity; }
  if(fValueCT) { delete [] fValueCT; }
  if(fMergedSlices) { delete fMergedSlices; }
  
}

GetDicomDataPath()

G4String DicomHandler::GetDicomDataPath ( )

static

Definition at line 83 of file DicomHandler.cc.

{
    // default is current directory
    G4String driverPath = ".";
    // check environment
    const char* path = std::getenv("DICOM_PATH");
 
    if(path)
    {
        // if is set in environment
        return G4String(path);
    }
    else
    {
        // if DICOM_USE_HEAD, look for data installation
#ifdef DICOM_USE_HEAD
        G4cerr << "Warning! DICOM was compiled with DICOM_USE_HEAD option but "
               << "the DICOM_PATH was not set!" << G4endl;
        G4String datadir = G4GetEnv<G4String>("G4ENSDFSTATEDATA", "");
        if(datadir.length() > 0)
        {
            auto _last = datadir.rfind("/");
            if(_last != std::string::npos)
                datadir.erase(_last);
            driverPath = datadir + "/DICOM1.1/DICOM_HEAD_TEST";
            G4int rc = setenv("DICOM_PATH", driverPath.c_str(), 0);
            G4cerr << "\t --> Using '" << driverPath << "'..." << G4endl;
            G4ConsumeParameters(rc);
        }
#endif
    }
    return driverPath;
}

GetDicomDataFile()

G4String DicomHandler::GetDicomDataFile ( )

static

Definition at line 119 of file DicomHandler.cc.

{
#if defined(DICOM_USE_HEAD) && defined(G4_DCMTK)
    return GetDicomDataPath() + "/Data.dat.new";
#else
    return GetDicomDataPath() + "/Data.dat";
#endif
}

GetValue()

template<class Type >

void DicomHandler::GetValue ( char *  _val, Type &  _rval  )

private

Definition at line 1204 of file DicomHandler.cc.

                                                     {
  
#if BYTE_ORDER == BIG_ENDIAN
  if(fLittleEndian) {      // little endian
#else // BYTE_ORDER == LITTLE_ENDIAN
    if(!fLittleEndian) {     // big endian
#endif
      const G4int SIZE = sizeof(_rval);
      char ctemp;
      for(G4int i = 0; i < SIZE/2; ++i) {
        ctemp = _val[i];
        _val[i] = _val[SIZE - 1 - i];
        _val[SIZE - 1 - i] = ctemp;
      }
    }
    _rval = *(Type *)_val;
  }

ReadCalibration()

void DicomHandler::ReadCalibration ( )

private

Definition at line 876 of file DicomHandler.cc.

{
 fNbrequali = 0;
// CT2Density.dat contains the calibration curve to convert CT (Hounsfield)
// number to physical density
 std::ifstream calibration(fCt2DensityFile.c_str());
 calibration >> fNbrequali; 
 fValueDensity = new G4double[fNbrequali];
 fValueCT = new G4double[fNbrequali];
 
 if(!calibration) {
   G4Exception("DicomHandler::ReadFile","DICOM001", FatalException,
               "@@@ No value to transform pixels in density!");
   } 
   else { // calibration was successfully opened
      for(G4int i = 0; i < fNbrequali; ++i) 
         { // Loop to store all the pts in CT2Density.dat
        calibration >> fValueCT[i] >> fValueDensity[i];
        }
  }
calibration.close();
fReadCalibration = true;
}

GetInformation()

void DicomHandler::GetInformation ( G4int &  tagDictionary, char *  data  )

private

Definition at line 361 of file DicomHandler.cc.

{
    if(tagDictionary == 0x00280010 ) { // Number of Rows
        GetValue(data, fRows);
        std::printf("[0x00280010] Rows -> %i\n",fRows);
 
    } else if(tagDictionary == 0x00280011 ) { // Number of fColumns
        GetValue(data, fColumns);
        std::printf("[0x00280011] Columns -> %i\n",fColumns);
 
    } else if(tagDictionary == 0x00280102 ) { // High bits  ( not used )
        short highBits;
        GetValue(data, highBits);
        std::printf("[0x00280102] High bits -> %i\n",highBits);
 
    } else if(tagDictionary == 0x00280100 ) { // Bits allocated
        GetValue(data, fBitAllocated);
        std::printf("[0x00280100] Bits allocated -> %i\n", fBitAllocated);
 
    } else if(tagDictionary == 0x00280101 ) { //  Bits stored ( not used )
        short bitStored;
        GetValue(data, bitStored);
        std::printf("[0x00280101] Bits stored -> %i\n",bitStored);
 
    } else if(tagDictionary == 0x00280106 ) { //  Min. pixel value
        GetValue(data, fMinPixelValue);
        std::printf("[0x00280106] Min. pixel value -> %i\n", fMinPixelValue);
 
    } else if(tagDictionary == 0x00280107 ) { //  Max. pixel value
        GetValue(data, fMaxPixelValue);
        std::printf("[0x00280107] Max. pixel value -> %i\n", fMaxPixelValue);
 
    } else if(tagDictionary == 0x00281053) { //  Rescale slope
        fRescaleSlope = atoi(data);
        std::printf("[0x00281053] Rescale Slope -> %d\n", fRescaleSlope);
 
    } else if(tagDictionary == 0x00281052 ) { // Rescalse intercept
        fRescaleIntercept = atoi(data);
        std::printf("[0x00281052] Rescale Intercept -> %d\n", 
                    fRescaleIntercept );
 
    } else if(tagDictionary == 0x00280103 ) {
        //  Pixel representation ( functions not design to read signed bits )
      fPixelRepresentation = atoi(data); // 0: unsigned  1: signed
        std::printf("[0x00280103] Pixel Representation -> %i\n",
                    fPixelRepresentation);
        if(fPixelRepresentation == 1 ) {
            std::printf("### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
            std::printf("DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
            std::printf("ERROR !!!!!! -> \n");
        }
 
    } else if(tagDictionary == 0x00080006 ) { //  Modality
        std::printf("[0x00080006] Modality -> %s\n", data);
 
    } else if(tagDictionary == 0x00080070 ) { //  Manufacturer
        std::printf("[0x00080070] Manufacturer -> %s\n", data);
 
    } else if(tagDictionary == 0x00080080 ) { //  Institution Name
        std::printf("[0x00080080] Institution Name -> %s\n", data);
 
    } else if(tagDictionary == 0x00080081 ) { //  Institution Address
        std::printf("[0x00080081] Institution Address -> %s\n", data);
 
    } else if(tagDictionary == 0x00081040 ) { //  Institution Department Name
        std::printf("[0x00081040] Institution Department Name -> %s\n", data);
 
    } else if(tagDictionary == 0x00081090 ) { //  Manufacturer's Model Name
        std::printf("[0x00081090] Manufacturer's Model Name -> %s\n", data);
 
    } else if(tagDictionary == 0x00181000 ) { //  Device Serial Number
        std::printf("[0x00181000] Device Serial Number -> %s\n", data);
 
    } else if(tagDictionary == 0x00080008 ) { //  Image type ( not used )
        std::printf("[0x00080008] Image Types -> %s\n", data);
 
    } else if(tagDictionary == 0x00283000 ) { //Modality LUT Sequence(not used)
        std::printf("[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
 
    } else if(tagDictionary == 0x00283002 ) { // LUT Descriptor ( not used )
        std::printf("[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
 
    } else if(tagDictionary == 0x00283003 ) { // LUT Explanation ( not used )
        std::printf("[0x00283003] LUT Explanation LO 1 -> %s\n", data);
 
    } else if(tagDictionary == 0x00283004 ) { // Modality LUT ( not used )
        std::printf("[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
 
    } else if(tagDictionary == 0x00283006 ) { // LUT Data ( not used )
        std::printf("[0x00283006] LUT Data US or SS -> %s\n", data);
 
    } else if(tagDictionary == 0x00283010 ) { // VOI LUT ( not used )
        std::printf("[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
 
    } else if(tagDictionary == 0x00280120 ) { // Pixel Padding Value (not used)
      std::printf("[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
 
    } else if(tagDictionary == 0x00280030 ) { // Pixel Spacing
        G4String datas(data);
        G4int iss = G4int(datas.find('\\'));
        fPixelSpacingX = atof( datas.substr(0,iss).c_str() );
        fPixelSpacingY = atof( datas.substr(iss+1,datas.length()).c_str() );
 
    } else if(tagDictionary == 0x00200037 ) { // Image Orientation ( not used )
        std::printf("[0x00200037] Image Orientation (Phantom) -> %s\n", data);
 
    } else if(tagDictionary == 0x00200032 ) { // Image Position ( not used )
        std::printf("[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
 
    } else if(tagDictionary == 0x00180050 ) { // Slice Thickness
      fSliceThickness = atof(data);
      std::printf("[0x00180050] Slice Thickness (mm) -> %f\n", fSliceThickness);
 
    } else if(tagDictionary == 0x00201041 ) { // Slice Location
      fSliceLocation = atof(data);
      std::printf("[0x00201041] Slice Location -> %f\n", fSliceLocation);
 
    } else if(tagDictionary == 0x00280004 ) { // Photometric Interpretation
      // ( not used )
        std::printf("[0x00280004] Photometric Interpretation -> %s\n", data);
 
    } else if(tagDictionary == 0x00020010) { // Endian
        if(strcmp(data, "1.2.840.10008.1.2") == 0)
            fImplicitEndian = true;
        else if(strncmp(data, "1.2.840.10008.1.2.2", 19) == 0)
            fLittleEndian = false;
        //else 1.2.840..10008.1.2.1 (explicit little endian)
 
        std::printf("[0x00020010] Endian -> %s\n", data);
    }
 
    // others
    else {
        //std::printf("[0x%x] -> %s\n", tagDictionary, data);
        ;
    }
 
}

Pixel2density()

G4float DicomHandler::Pixel2density ( G4int  pixel )

private

Definition at line 929 of file DicomHandler.cc.

{
  if(!fReadCalibration) { ReadCalibration(); }
  
  G4float density = -1.;
  G4double deltaCT = 0;
  G4double deltaDensity = 0;
  
  
  for(G4int j = 1; j < fNbrequali; ++j) {
    if( pixel >= fValueCT[j-1] && pixel < fValueCT[j]) {
      
      deltaCT = fValueCT[j] - fValueCT[j-1];
      deltaDensity = fValueDensity[j] - fValueDensity[j-1];
      
      // interpolating linearly
      density = G4float(fValueDensity[j]
                        -((fValueCT[j] - pixel)*deltaDensity/deltaCT ));
      break;
    }
  }
  
  if(density < 0.) {
    std::printf("@@@ Error density = %f && Pixel = %i \
      (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
                deltaDensity/deltaCT);
  }
  
  return density;
}

ReadMaterialIndices()

void DicomHandler::ReadMaterialIndices ( std::ifstream &  finData )

private

Definition at line 640 of file DicomHandler.cc.

{
  unsigned int nMate;
  G4String mateName;
  G4float densityMax;
  finData >> nMate;
  if( finData.eof() ) return;
  
  G4cout << " ReadMaterialIndices " << nMate << G4endl;
  for( unsigned int ii = 0; ii < nMate; ++ii )
  {
    finData >> mateName >> densityMax;
    fMaterialIndices[densityMax] = mateName;
    //    G4cout << ii << " ReadMaterialIndices " << mateName << " " 
    //     << densityMax << G4endl;
  }
  
}

GetMaterialIndex()

unsigned int DicomHandler::GetMaterialIndex ( G4float  density )

private

Definition at line 661 of file DicomHandler.cc.

{
  std::map<G4float,G4String>::const_reverse_iterator ite;
  std::size_t ii = fMaterialIndices.size();
 
  for( ite = fMaterialIndices.crbegin(); ite != fMaterialIndices.crend();
       ++ite, ii-- )
  {
    if( density >= (*ite).first )  { break; }
  }
  
  if(ii == fMaterialIndices.size())  { ii = fMaterialIndices.size()-1; }
  
  return unsigned(ii);
}

StoreData() [1/2]

void DicomHandler::StoreData ( std::ofstream &  foutG4DCM )

private

Definition at line 556 of file DicomHandler.cc.

{
  G4int mean;
  G4double density;
  G4bool overflow = false;
  
  //----- Print indices of material
  if(fCompression == 1) { // no fCompression: each pixel has a density value)
    for( G4int ww = 0; ww < fRows; ++ww) {
      for( G4int xx = 0; xx < fColumns; ++xx) {
        mean = fTab[ww][xx];
        density = Pixel2density(mean);
        foutG4DCM << GetMaterialIndex( G4float(density) ) << " ";
      }
      foutG4DCM << G4endl;
    }
    
  } else {
    // density value is the average of a square region of
    // fCompression*fCompression pixels
    for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
      for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
        overflow = false;
        mean = 0;
        for(G4int sumx = 0; sumx < fCompression; ++sumx) {
          for(G4int sumy = 0; sumy < fCompression; ++sumy) {
            if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
            mean += fTab[ww+sumy][xx+sumx];
          }
          if(overflow) break;
        }
        mean /= fCompression*fCompression;
        
        if(!overflow) {
          density = Pixel2density(mean);
          foutG4DCM << GetMaterialIndex( G4float(density) ) << " ";
        }
      }
      foutG4DCM << G4endl;
    }
    
  }
  
  //----- Print densities
  if(fCompression == 1) { // no fCompression: each pixel has a density value)
    for( G4int ww = 0; ww < fRows; ww++) {
      for( G4int xx = 0; xx < fColumns; xx++) {
        mean = fTab[ww][xx];
        density = Pixel2density(mean);
        foutG4DCM << density << " ";
        if( xx%8 == 3 ) foutG4DCM << G4endl; // just for nicer reading
      }
    }
    
  } else {
    // density value is the average of a square region of
    // fCompression*fCompression pixels
    for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
      for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
        overflow = false;
        mean = 0;
        for(G4int sumx = 0; sumx < fCompression; ++sumx) {
          for(G4int sumy = 0; sumy < fCompression; ++sumy) {
            if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
            mean += fTab[ww+sumy][xx+sumx];
          }
          if(overflow) break;
        }
        mean /= fCompression*fCompression;
        
        if(!overflow) {
          density = Pixel2density(mean);
          foutG4DCM << density  << " ";
          if( xx/fCompression%8 == 3 ) foutG4DCM << G4endl; // just for nicer
          // reading
        }
      }
    }
    
  }
  
}

StoreData() [2/2]

void DicomHandler::StoreData ( DicomPhantomZSliceHeader *  dcmPZSH )

private

Definition at line 502 of file DicomHandler.cc.

{
    G4int mean;
    G4double density;
    G4bool overflow = false;
 
    if(!dcmPZSH) { return; }
 
    dcmPZSH->SetSliceLocation(fSliceLocation);
 
    //----- Print indices of material
    if(fCompression == 1) { // no fCompression: each pixel has a density value)
        for( G4int ww = 0; ww < fRows; ++ww) {
            dcmPZSH->AddRow();
            for( G4int xx = 0; xx < fColumns; ++xx) {
                mean = fTab[ww][xx];
                density = Pixel2density(mean);
                dcmPZSH->AddValue(density);
                dcmPZSH->AddMateID(GetMaterialIndex(G4float(density)));
            }
        }
 
    } else {
        // density value is the average of a square region of
        // fCompression*fCompression pixels
      for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
        dcmPZSH->AddRow();
        for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
          overflow = false;
          mean = 0;
          for(G4int sumx = 0; sumx < fCompression; ++sumx) {
            for(G4int sumy = 0; sumy < fCompression; ++sumy) {
              if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
              mean += fTab[ww+sumy][xx+sumx];
            }
            if(overflow) break;
          }
          mean /= fCompression*fCompression;
          
          if(!overflow) {
            density = Pixel2density(mean);
            dcmPZSH->AddValue(density);
            dcmPZSH->AddMateID(GetMaterialIndex(G4float(density)));
          }
        }
      }
    }
    
    dcmPZSH->FlipData();
}

read_defined_nested()

G4int DicomHandler::read_defined_nested ( FILE *  nested, G4int  SQ_Length  )

private

Definition at line 1098 of file DicomHandler.cc.

{
  //      VARIABLES
  unsigned short item_GroupNumber;
  unsigned short item_ElementNumber;
  G4int item_Length;
  G4int items_array_length=0;
  char * buffer= new char[LINEBUFFSIZE];
  size_t rflag = 0;
  
  while(items_array_length < SQ_Length)
    {
      rflag = std::fread(buffer, 2, 1, nested);
      GetValue(buffer, item_GroupNumber);
      
      rflag = std::fread(buffer, 2, 1, nested);
      GetValue(buffer, item_ElementNumber);
      
      rflag = std::fread(buffer, 4, 1, nested);
      GetValue(buffer, item_Length);
      
      rflag = std::fread(buffer, item_Length, 1, nested);
      
      items_array_length= items_array_length+8+item_Length;
    }
  
  delete [] buffer;
  
  if( SQ_Length>items_array_length )
    return 0;
  else
    return 1;
  if (rflag) return 1;
}

read_undefined_nested()

void DicomHandler::read_undefined_nested ( FILE *  nested )

private

Definition at line 1135 of file DicomHandler.cc.

{
  //      VARIABLES
  unsigned short item_GroupNumber;
  unsigned short item_ElementNumber;
  unsigned int item_Length;
  char * buffer= new char[LINEBUFFSIZE];
  size_t rflag = 0;
  
    do
      {
        rflag = std::fread(buffer, 2, 1, nested);
        GetValue(buffer, item_GroupNumber);
        
        rflag = std::fread(buffer, 2, 1, nested);
        GetValue(buffer, item_ElementNumber);
        
        rflag = std::fread(buffer, 4, 1, nested);
        GetValue(buffer, item_Length);
        
        if(item_Length!=0xffffffff)
          rflag = std::fread(buffer, item_Length, 1, nested);
        else
          read_undefined_item(nested);
        
        
      } while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD 
              || item_Length!=0);
    
    delete [] buffer;
    if (rflag) return;
}

read_undefined_item()

void DicomHandler::read_undefined_item ( FILE *  nested )

private

Definition at line 1170 of file DicomHandler.cc.

{
  //      VARIABLES
  unsigned short item_GroupNumber;
  unsigned short item_ElementNumber;
  G4int item_Length; size_t rflag = 0;
  char *buffer= new char[LINEBUFFSIZE];
  
  do
    {
      rflag = std::fread(buffer, 2, 1, nested);
      GetValue(buffer, item_GroupNumber);
      
      rflag = std::fread(buffer, 2, 1, nested);
      GetValue(buffer, item_ElementNumber);
      
      rflag = std::fread(buffer, 4, 1, nested);
      GetValue(buffer, item_Length);
      
      
      if(item_Length!=0)
        rflag = std::fread(buffer,item_Length,1,nested);
      
    }
  while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D 
        || item_Length!=0);
  
  delete [] buffer;
  if (rflag) return;
}

Member Data Documentation

fInstance

DicomHandler * DicomHandler::fInstance = 0

staticprivate

DicomHandler.cc :

• Handling of DICM images
  • Reading headers and pixels
• Transforming pixel to density and creating *.g4dcm files

Definition at line 97 of file DicomHandler.hh.

DATABUFFSIZE

const G4int DicomHandler::DATABUFFSIZE

private

Definition at line 99 of file DicomHandler.hh.

LINEBUFFSIZE

const G4int DicomHandler::LINEBUFFSIZE

private

Definition at line 100 of file DicomHandler.hh.

FILENAMESIZE

const G4int DicomHandler::FILENAMESIZE

private

Definition at line 101 of file DicomHandler.hh.

fCompression

short DicomHandler::fCompression

private

Definition at line 114 of file DicomHandler.hh.

fNFiles

G4int DicomHandler::fNFiles

private

Definition at line 115 of file DicomHandler.hh.

fRows

short DicomHandler::fRows

private

Definition at line 116 of file DicomHandler.hh.

fColumns

short DicomHandler::fColumns

private

Definition at line 117 of file DicomHandler.hh.

fBitAllocated

short DicomHandler::fBitAllocated

private

Definition at line 118 of file DicomHandler.hh.

fMaxPixelValue

G4int DicomHandler::fMaxPixelValue

private

Definition at line 119 of file DicomHandler.hh.

fMinPixelValue

G4int DicomHandler::fMinPixelValue

private

Definition at line 119 of file DicomHandler.hh.

fPixelSpacingX

G4double DicomHandler::fPixelSpacingX

private

Definition at line 121 of file DicomHandler.hh.

fPixelSpacingY

G4double DicomHandler::fPixelSpacingY

private

Definition at line 121 of file DicomHandler.hh.

fSliceThickness

G4double DicomHandler::fSliceThickness

private

Definition at line 122 of file DicomHandler.hh.

fSliceLocation

G4double DicomHandler::fSliceLocation

private

Definition at line 123 of file DicomHandler.hh.

fRescaleIntercept

G4int DicomHandler::fRescaleIntercept

private

Definition at line 125 of file DicomHandler.hh.

fRescaleSlope

G4int DicomHandler::fRescaleSlope

private

Definition at line 125 of file DicomHandler.hh.

fLittleEndian

G4bool DicomHandler::fLittleEndian

private

Definition at line 127 of file DicomHandler.hh.

fImplicitEndian

G4bool DicomHandler::fImplicitEndian

private

Definition at line 127 of file DicomHandler.hh.

fPixelRepresentation

short DicomHandler::fPixelRepresentation

private

Definition at line 128 of file DicomHandler.hh.

fTab

G4int** DicomHandler::fTab

private

Definition at line 130 of file DicomHandler.hh.

fMaterialIndices

std::map<G4float,G4String> DicomHandler::fMaterialIndices

private

Definition at line 131 of file DicomHandler.hh.

fNbrequali

G4int DicomHandler::fNbrequali

private

Definition at line 133 of file DicomHandler.hh.

fValueDensity

G4double* DicomHandler::fValueDensity

private

Definition at line 134 of file DicomHandler.hh.

fValueCT

G4double* DicomHandler::fValueCT

private

Definition at line 135 of file DicomHandler.hh.

fReadCalibration

G4bool DicomHandler::fReadCalibration

private

Definition at line 136 of file DicomHandler.hh.

fMergedSlices

DicomPhantomZSliceMerged* DicomHandler::fMergedSlices

private

Definition at line 137 of file DicomHandler.hh.

fDriverFile

G4String DicomHandler::fDriverFile

private

Definition at line 139 of file DicomHandler.hh.

fCt2DensityFile

G4String DicomHandler::fCt2DensityFile

private

Definition at line 140 of file DicomHandler.hh.

---

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

• .doxygen/Doxymodules_medical.h
• extended/medical/DICOM/include/DicomHandler.hh
• extended/medical/DICOM/src/DicomHandler.cc