DicomHandler Class Reference

#include <DicomHandler.hh>

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

Static Public Member Functions
static DicomHandler *	Instance ()

Detailed Description

Definition at line 71 of file DicomHandler.hh.

Constructor & Destructor Documentation

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), nbrequali(0),
    valuedensity(NULL),valueCT(NULL),readCalibration(false),
    mergedSlices(NULL),driverFile("Data.dat"),ct2densityFile("CT2Density.dat")
{
    mergedSlices = new DicomPhantomZSliceMerged;
}

DicomHandler::~DicomHandler

(

)

Definition at line 97 of file DicomHandler.cc.

{

}

Member Function Documentation

void DicomHandler::CheckFileFormat

(

)

Definition at line 871 of file DicomHandler.cc.

References DicomPhantomZSliceMerged::CheckSlices(), G4String::data(), G4cout, G4endl, and ReadFile().

Referenced by main().

{
    std::ifstream checkData(driverFile.c_str());
    char * oneLine = new char[128];

    if(!(checkData.is_open())) { //Check existance of Data.dat

        G4cout << "\nDicomG4 needs Data.dat (or another driver file specified in command line)"
        << ":\n\tFirst line: number of image pixel for a "
        << "voxel (G4Box)\n\tSecond line: number of images (CT slices) to "
        << "read\n\tEach following line contains the name of a Dicom image except "
        << "for the .dcm extension\n";
        exit(0);
    }

    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(driverFile.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);
            std::cout << "check 1: " << name << std::endl;
            //  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
        mergedSlices->CheckSlices();

        delete [] fCompressionc;
        delete [] maxc;
        delete [] name;
        delete [] inputFile;
        if (rflag) return;

    }

    if(valuedensity) { delete [] valuedensity; }
    if(valueCT) { delete [] valueCT; }
    if(mergedSlices) { delete mergedSlices; }


}

DicomHandler * DicomHandler::Instance ( void  )

static

Definition at line 72 of file DicomHandler.cc.

{
    return fgInstance;
}

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

Definition at line 624 of file DicomHandler.cc.

References density.

Referenced by ReadFile().

{
    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;
}

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

Definition at line 104 of file DicomHandler.cc.

References DicomPhantomZSliceHeader::AddMaterial(), DicomPhantomZSliceMerged::AddZSlice(), buffer, G4cerr, G4cout, G4endl, ReadData(), DicomPhantomZSliceHeader::SetMaxX(), DicomPhantomZSliceHeader::SetMaxY(), DicomPhantomZSliceHeader::SetMaxZ(), DicomPhantomZSliceHeader::SetMinX(), DicomPhantomZSliceHeader::SetMinY(), DicomPhantomZSliceHeader::SetMinZ(), DicomPhantomZSliceHeader::SetNoVoxelX(), DicomPhantomZSliceHeader::SetNoVoxelY(), and DicomPhantomZSliceHeader::SetNoVoxelZ().

Referenced by CheckFileFormat().

{
    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){
                        G4cerr << "Function read_defined_nested() failed!" << G4endl;
                        exit(-10);               }
                }
            } 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) {
                    G4cerr << "[DicomHandler] fseek failed" << G4endl;
                    exit(-10);}

                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.);

    //=====================================================================
    // This is depreciated --> Handled by DicomPhantomZSliceHeader
    /*
     // Creating files to store information
     std::ofstream foutG4DCM;
     foutG4DCM.open(fnameG4DCM);
     G4cout << "### Writing of " << fnameG4DCM << " ### " << G4endl;

     foutG4DCM << fMaterialIndices.size() << G4endl;
     //--- Write materials
     unsigned int ii = 0;
     for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++, ii++ ){
     foutG4DCM << ii << " " << (*ite).second << G4endl;
     }
     //--- Write number of voxels (assume only one voxel in Z)
     foutG4DCM << fColumns/fCompression << " " << fRows/fCompression << " 1 " << G4endl;
     //--- Write minimum and maximum extensions
     foutG4DCM << -fPixelSpacingX*fColumns/2. << " " << fPixelSpacingX*fColumns/2. << G4endl;
     foutG4DCM << -fPixelSpacingY*fRows/2. << " " << fPixelSpacingY*fRows/2. <<
     G4endl;
     foutG4DCM << fSliceLocation-fSliceThickness/2. << " "
     << fSliceLocation+fSliceThickness/2. << G4endl;
     //    foutG4DCM << fCompression << G4endl;
     */
    //=====================================================================

    ReadData( dicom, filename2 );

    // DEPRECIATED
    //StoreData( foutG4DCM );
    //foutG4DCM.close();

    StoreData( zslice );

    // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency
    //zslice->DumpToFile();

    mergedSlices->AddZSlice(zslice);

    //
    delete [] buffer;
    delete [] data;

    if (rflag) return returnvalue;
    return returnvalue;
}

---

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

• DicomHandler.hh
• DicomHandler.cc