DicomHandler Class Reference

#include <DicomHandler.hh>

Collaboration diagram for DicomHandler:

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

Static Public Member Functions
static DicomHandler *	Instance ()

Private Member Functions
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 int	DATABUFFSIZE
const int	LINEBUFFSIZE
const int	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
bool	fReadCalibration
DicomPhantomZSliceMerged *	fMergedSlices
G4String	fDriverFile
G4String	fCt2DensityFile

Static Private Attributes
static DicomHandler *	fInstance = 0

Detailed Description

Definition at line 71 of file DicomHandler.hh.

Constructor & Destructor Documentation

DicomHandler()

DicomHandler::DicomHandler

(

)

Definition at line 79 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;
}

~DicomHandler()

DicomHandler::~DicomHandler

(

)

Definition at line 96 of file DicomHandler.cc.

{

}

Member Function Documentation

CheckFileFormat()

void DicomHandler::CheckFileFormat

(

)

Definition at line 863 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 * name = new char[FILENAMESIZE];
    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;
    
    for( G4int i = 1; i <= fNFiles; i++ ) { // Begin loop on filenames
      
      rflag = std::fscanf(lecturePref,"%s",inputFile);
      std::sprintf(name,"%s.dcm",inputFile);
      G4cout << "check 1: " << name << G4endl;
      //  Open input file and give it to gestion_dicom :
      std::printf("### Opening %s and reading :\n",name);
      dicom = std::fopen(name,"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);
    }
    rflag = std::fclose(lecturePref);
    
    // Checks the spacing is correct. Dumps to files
    fMergedSlices->CheckSlices();
    
    delete [] fCompressionc;
    delete [] maxc;
    delete [] name;
    delete [] inputFile;
    if (rflag) return;
    
  }
  
  if(fValueDensity) { delete [] fValueDensity; }
  if(fValueCT) { delete [] fValueCT; }
  if(fMergedSlices) { delete fMergedSlices; }
  
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetInformation()

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

private

Definition at line 293 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);
        int iss = datas.find('\\');
        fPixelSpacingX = atof( datas.substr(0,iss).c_str() );
        fPixelSpacingY = atof( datas.substr(iss+2,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);
        ;
    }

}

Here is the call graph for this function:

Here is the caller graph for this function:

GetMaterialIndex()

unsigned int DicomHandler::GetMaterialIndex ( G4float  density )

private

Definition at line 592 of file DicomHandler.cc.

{
  std::map<G4float,G4String>::reverse_iterator ite;
  G4int ii = fMaterialIndices.size();
  for( ite = fMaterialIndices.rbegin(); ite != fMaterialIndices.rend();
       ite++, ii-- ) {
    if( density >= (*ite).first ) {
      break;
    }
  }
  //-  G4cout << " GetMaterialIndex " << density << " = " << ii << G4endl;
  
  if(static_cast<unsigned int>(ii) == fMaterialIndices.size())
    { ii = fMaterialIndices.size()-1; }
  
  return  ii;
  
}

Here is the caller graph for this function:

GetValue()

template<class Type >

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

private

Definition at line 1075 of file DicomHandler.cc.

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

Here is the caller graph for this function:

Instance()

DicomHandler * DicomHandler::Instance ( void  )

static

Definition at line 72 of file DicomHandler.cc.

{
    return fInstance;
}

Pixel2density()

G4float DicomHandler::Pixel2density ( G4int  pixel )

private

Definition at line 831 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 = 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

read_defined_nested()

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

private

Definition at line 969 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

read_undefined_item()

void DicomHandler::read_undefined_item ( FILE *  nested )

private

Definition at line 1041 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

read_undefined_nested()

void DicomHandler::read_undefined_nested ( FILE *  nested )

private

Definition at line 1006 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadCalibration()

void DicomHandler::ReadCalibration ( )

private

Definition at line 800 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadData()

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

Definition at line 613 of file DicomHandler.cc.

{
  G4int returnvalue = 0; size_t rflag = 0;
  
  //  READING THE PIXELS :
  G4int w = 0;
  
  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++) {
        w++;
        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++) {
        w++;
        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
  char * nameProcessed = new char[FILENAMESIZE];
  FILE* fileOut;
  
  std::sprintf(nameProcessed,"%s.g4dcmb",filename2);
  fileOut = std::fopen(nameProcessed,"w+b");
  std::printf("### Writing of %s ###\n",nameProcessed);
  
  unsigned int nMate = fMaterialIndices.size();
  rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
  //--- Write materials
  std::map<G4float,G4String>::const_iterator ite;
  for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++ ){
    G4String mateName = (*ite).second;
    for( G4int 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 = -fPixelSpacingX*fColumns/2.;
  G4float pixelLocationXP = fPixelSpacingX*fColumns/2.;
  G4float pixelLocationYM = -fPixelSpacingY*fRows/2.;
  G4float pixelLocationYP = fPixelSpacingY*fRows/2.;
  G4float fSliceLocationZM = fSliceLocation-fSliceThickness/2.;
  G4float fSliceLocationZP = 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(int sumx = 0; sumx < compSize; sumx++) {
          for(int 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(int sumx = 0; sumx < compSize; sumx++) {
          for(int 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadFile()

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

Definition at line 103 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) {
        // Folling 2 fread's are modifications to original DICOM example 
        //(Jonathan Madsen)
        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);

    std::map<G4float,G4String>::const_iterator ite;
    for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ++ite){
      zslice->AddMaterial(ite->second);
    }
    
    zslice->SetNoVoxelX(fColumns/fCompression);
    zslice->SetNoVoxelY(fRows/fCompression);
    zslice->SetNoVoxelZ(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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadMaterialIndices()

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

private

Definition at line 572 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;
  }
  
}

Here is the caller graph for this function:

StoreData() [1/2]

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

private

Definition at line 488 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( 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(int sumx = 0; sumx < fCompression; sumx++) {
          for(int 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( 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(int sumx = 0; sumx < fCompression; sumx++) {
          for(int 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
        }
      }
    }
    
  }
  
}

Here is the call graph for this function:

Here is the caller graph for this function:

StoreData() [2/2]

void DicomHandler::StoreData ( DicomPhantomZSliceHeader *  dcmPZSH )

private

Definition at line 434 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(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(int sumx = 0; sumx < fCompression; sumx++) {
            for(int 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(density));
          }
        }
      }
    }
    
    dcmPZSH->FlipData();
}

Here is the call graph for this function:

Member Data Documentation

DATABUFFSIZE

const int DicomHandler::DATABUFFSIZE

private

Definition at line 97 of file DicomHandler.hh.

fBitAllocated

short DicomHandler::fBitAllocated

private

Definition at line 116 of file DicomHandler.hh.

fColumns

short DicomHandler::fColumns

private

Definition at line 115 of file DicomHandler.hh.

fCompression

short DicomHandler::fCompression

private

Definition at line 112 of file DicomHandler.hh.

fCt2DensityFile

G4String DicomHandler::fCt2DensityFile

private

Definition at line 138 of file DicomHandler.hh.

fDriverFile

G4String DicomHandler::fDriverFile

private

Definition at line 137 of file DicomHandler.hh.

FILENAMESIZE

const int DicomHandler::FILENAMESIZE

private

Definition at line 99 of file DicomHandler.hh.

fImplicitEndian

G4bool DicomHandler::fImplicitEndian

private

Definition at line 125 of file DicomHandler.hh.

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 95 of file DicomHandler.hh.

fLittleEndian

G4bool DicomHandler::fLittleEndian

private

Definition at line 125 of file DicomHandler.hh.

fMaterialIndices

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

private

Definition at line 129 of file DicomHandler.hh.

fMaxPixelValue

G4int DicomHandler::fMaxPixelValue

private

Definition at line 117 of file DicomHandler.hh.

fMergedSlices

DicomPhantomZSliceMerged* DicomHandler::fMergedSlices

private

Definition at line 135 of file DicomHandler.hh.

fMinPixelValue

G4int DicomHandler::fMinPixelValue

private

Definition at line 117 of file DicomHandler.hh.

fNbrequali

G4int DicomHandler::fNbrequali

private

Definition at line 131 of file DicomHandler.hh.

fNFiles

G4int DicomHandler::fNFiles

private

Definition at line 113 of file DicomHandler.hh.

fPixelRepresentation

short DicomHandler::fPixelRepresentation

private

Definition at line 126 of file DicomHandler.hh.

fPixelSpacingX

G4double DicomHandler::fPixelSpacingX

private

Definition at line 119 of file DicomHandler.hh.

fPixelSpacingY

G4double DicomHandler::fPixelSpacingY

private

Definition at line 119 of file DicomHandler.hh.

fReadCalibration

bool DicomHandler::fReadCalibration

private

Definition at line 134 of file DicomHandler.hh.

fRescaleIntercept

G4int DicomHandler::fRescaleIntercept

private

Definition at line 123 of file DicomHandler.hh.

fRescaleSlope

G4int DicomHandler::fRescaleSlope

private

Definition at line 123 of file DicomHandler.hh.

fRows

short DicomHandler::fRows

private

Definition at line 114 of file DicomHandler.hh.

fSliceLocation

G4double DicomHandler::fSliceLocation

private

Definition at line 121 of file DicomHandler.hh.

fSliceThickness

G4double DicomHandler::fSliceThickness

private

Definition at line 120 of file DicomHandler.hh.

fTab

G4int** DicomHandler::fTab

private

Definition at line 128 of file DicomHandler.hh.

fValueCT

G4double* DicomHandler::fValueCT

private

Definition at line 133 of file DicomHandler.hh.

fValueDensity

G4double* DicomHandler::fValueDensity

private

Definition at line 132 of file DicomHandler.hh.

LINEBUFFSIZE

const int DicomHandler::LINEBUFFSIZE

private

Definition at line 98 of file DicomHandler.hh.

---

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

• DicomHandler.hh
• DicomHandler.cc