Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/DICOM/src/DicomHandler.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/DICOM/src/DicomHandler.cc (Version 11.3.0) and /examples/extended/medical/DICOM/src/DicomHandler.cc (Version 10.2.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: DicomHandler.cc 92820 2015-09-17 15:22:14Z gcosmo $
 26 //                                                 27 //
 27 /// \file medical/DICOM/src/DicomHandler.cc        28 /// \file medical/DICOM/src/DicomHandler.cc
 28 /// \brief Implementation of the DicomHandler      29 /// \brief Implementation of the DicomHandler class
 29 //                                                 30 //
 30 // The code was written by :                       31 // The code was written by :
 31 //      *Louis Archambault louis.archambault@p     32 //      *Louis Archambault louis.archambault@phy.ulaval.ca,
 32 //      *Luc Beaulieu beaulieu@phy.ulaval.ca       33 //      *Luc Beaulieu beaulieu@phy.ulaval.ca
 33 //      +Vincent Hubert-Tremblay at tigre.2@sy     34 //      +Vincent Hubert-Tremblay at tigre.2@sympatico.ca
 34 //                                                 35 //
 35 //                                                 36 //
 36 // *Centre Hospitalier Universitaire de Quebec     37 // *Centre Hospitalier Universitaire de Quebec (CHUQ),
 37 // Hotel-Dieu de Quebec, departement de Radio-     38 // Hotel-Dieu de Quebec, departement de Radio-oncologie
 38 // 11 cote du palais. Quebec, QC, Canada, G1R      39 // 11 cote du palais. Quebec, QC, Canada, G1R 2J6
 39 // tel (418) 525-4444 #6720                        40 // tel (418) 525-4444 #6720
 40 // fax (418) 691 5268                              41 // fax (418) 691 5268
 41 //                                                 42 //
 42 // + University Laval, Quebec (QC) Canada          43 // + University Laval, Quebec (QC) Canada
 43 //********************************************     44 //*******************************************************
 44 //                                                 45 //
 45 //********************************************     46 //*******************************************************
 46 //                                                 47 //
 47 /// DicomHandler.cc :                              48 /// DicomHandler.cc :
 48 ///        - Handling of DICM images               49 ///        - Handling of DICM images
 49 ///         - Reading headers and pixels           50 ///         - Reading headers and pixels
 50 ///        - Transforming pixel to density and     51 ///        - Transforming pixel to density and creating *.g4dcm
 51 ///          files                                 52 ///          files
 52 //********************************************     53 //*******************************************************
 53                                                    54 
 54 #include "DicomHandler.hh"                         55 #include "DicomHandler.hh"
 55                                                << 
 56 #include "DicomPhantomZSliceHeader.hh"         << 
 57 #include "DicomPhantomZSliceMerged.hh"         << 
 58                                                << 
 59 #include "G4ios.hh"                            << 
 60 #include "globals.hh"                              56 #include "globals.hh"
                                                   >>  57 #include "G4ios.hh"
                                                   >>  58 #include <fstream>
 61                                                    59 
 62 #include <cctype>                                  60 #include <cctype>
 63 #include <cstring>                                 61 #include <cstring>
 64 #include <fstream>                             <<  62 
                                                   >>  63 #include "DicomPhantomZSliceHeader.hh"
                                                   >>  64 #include "DicomPhantomZSliceMerged.hh"
 65                                                    65 
 66 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67                                                    67 
 68 DicomHandler* DicomHandler::fInstance = 0;         68 DicomHandler* DicomHandler::fInstance = 0;
 69                                                    69 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 71                                                    71 
 72 DicomHandler* DicomHandler::Instance()             72 DicomHandler* DicomHandler::Instance()
 73 {                                                  73 {
 74   if (fInstance == 0) {                        <<  74     return fInstance;
 75     static DicomHandler dicomhandler;          << 
 76     fInstance = &dicomhandler;                 << 
 77   }                                            << 
 78   return fInstance;                            << 
 79 }                                                  75 }
 80                                                    76 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 82                                                    78 
 83 G4String DicomHandler::GetDicomDataPath()      <<  79 DicomHandler::DicomHandler()
                                                   >>  80 :   DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
                                                   >>  81     fCompression(0), fNFiles(0), fRows(0), fColumns(0),
                                                   >>  82     fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
                                                   >>  83     fPixelSpacingX(0.), fPixelSpacingY(0.),
                                                   >>  84     fSliceThickness(0.), fSliceLocation(0.),
                                                   >>  85     fRescaleIntercept(0), fRescaleSlope(0),
                                                   >>  86     fLittleEndian(true), fImplicitEndian(false),
                                                   >>  87     fPixelRepresentation(0), fNbrequali(0),
                                                   >>  88     fValueDensity(NULL),fValueCT(NULL),fReadCalibration(false),
                                                   >>  89     fMergedSlices(NULL),fDriverFile("Data.dat"),fCt2DensityFile("CT2Density.dat")
 84 {                                                  90 {
 85   // default is current directory              <<  91     fMergedSlices = new DicomPhantomZSliceMerged;
 86   G4String driverPath = ".";                   << 
 87   // check environment                         << 
 88   const char* path = std::getenv("DICOM_PATH") << 
 89                                                << 
 90   if (path) {                                  << 
 91     // if is set in environment                << 
 92     return G4String(path);                     << 
 93   }                                            << 
 94   else {                                       << 
 95     // if DICOM_USE_HEAD, look for data instal << 
 96 #ifdef DICOM_USE_HEAD                          << 
 97     G4cerr << "Warning! DICOM was compiled wit << 
 98            << "the DICOM_PATH was not set!" << << 
 99     G4String datadir = G4GetEnv<G4String>("G4E << 
100     if (datadir.length() > 0) {                << 
101       auto _last = datadir.rfind("/");         << 
102       if (_last != std::string::npos) datadir. << 
103       driverPath = datadir + "/DICOM1.1/DICOM_ << 
104       G4int rc = setenv("DICOM_PATH", driverPa << 
105       G4cerr << "\t --> Using '" << driverPath << 
106       G4ConsumeParameters(rc);                 << 
107     }                                          << 
108 #endif                                         << 
109   }                                            << 
110   return driverPath;                           << 
111 }                                                  92 }
112                                                    93 
113 //....oooOO0OOooo........oooOO0OOooo........oo     94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114                                                    95 
115 G4String DicomHandler::GetDicomDataFile()      <<  96 DicomHandler::~DicomHandler()
116 {                                                  97 {
117 #if defined(DICOM_USE_HEAD) && defined(G4_DCMT << 
118   return GetDicomDataPath() + "/Data.dat.new"; << 
119 #else                                          << 
120   return GetDicomDataPath() + "/Data.dat";     << 
121 #endif                                         << 
122 }                                              << 
123                                                << 
124 //....oooOO0OOooo........oooOO0OOooo........oo << 
125                                                    98 
126 #ifdef DICOM_USE_HEAD                          << 
127 DicomHandler::DicomHandler()                   << 
128   : DATABUFFSIZE(8192),                        << 
129     LINEBUFFSIZE(5020),                        << 
130     FILENAMESIZE(512),                         << 
131     fCompression(0),                           << 
132     fNFiles(0),                                << 
133     fRows(0),                                  << 
134     fColumns(0),                               << 
135     fBitAllocated(0),                          << 
136     fMaxPixelValue(0),                         << 
137     fMinPixelValue(0),                         << 
138     fPixelSpacingX(0.),                        << 
139     fPixelSpacingY(0.),                        << 
140     fSliceThickness(0.),                       << 
141     fSliceLocation(0.),                        << 
142     fRescaleIntercept(0),                      << 
143     fRescaleSlope(0),                          << 
144     fLittleEndian(true),                       << 
145     fImplicitEndian(false),                    << 
146     fPixelRepresentation(0),                   << 
147     fNbrequali(0),                             << 
148     fValueDensity(NULL),                       << 
149     fValueCT(NULL),                            << 
150     fReadCalibration(false),                   << 
151     fMergedSlices(NULL),                       << 
152     fCt2DensityFile("null.dat")                << 
153 {                                              << 
154   fMergedSlices = new DicomPhantomZSliceMerged << 
155   fDriverFile = GetDicomDataFile();            << 
156   G4cout << "Reading the DICOM_HEAD project "  << 
157 }                                                  99 }
158 #else                                          << 
159 DicomHandler::DicomHandler()                   << 
160   : DATABUFFSIZE(8192),                        << 
161     LINEBUFFSIZE(5020),                        << 
162     FILENAMESIZE(512),                         << 
163     fCompression(0),                           << 
164     fNFiles(0),                                << 
165     fRows(0),                                  << 
166     fColumns(0),                               << 
167     fBitAllocated(0),                          << 
168     fMaxPixelValue(0),                         << 
169     fMinPixelValue(0),                         << 
170     fPixelSpacingX(0.),                        << 
171     fPixelSpacingY(0.),                        << 
172     fSliceThickness(0.),                       << 
173     fSliceLocation(0.),                        << 
174     fRescaleIntercept(0),                      << 
175     fRescaleSlope(0),                          << 
176     fLittleEndian(true),                       << 
177     fImplicitEndian(false),                    << 
178     fPixelRepresentation(0),                   << 
179     fNbrequali(0),                             << 
180     fValueDensity(NULL),                       << 
181     fValueCT(NULL),                            << 
182     fReadCalibration(false),                   << 
183     fMergedSlices(NULL),                       << 
184     fDriverFile("Data.dat"),                   << 
185     fCt2DensityFile("CT2Density.dat")          << 
186 {                                              << 
187   fMergedSlices = new DicomPhantomZSliceMerged << 
188 }                                              << 
189 #endif                                         << 
190 //....oooOO0OOooo........oooOO0OOooo........oo << 
191                                                << 
192 DicomHandler::~DicomHandler() {}               << 
193                                                   100 
194 //....oooOO0OOooo........oooOO0OOooo........oo    101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
195                                                   102 
196 G4int DicomHandler::ReadFile(FILE* dicom, char    103 G4int DicomHandler::ReadFile(FILE* dicom, char* filename2)
197 {                                                 104 {
198   G4cout << " ReadFile " << filename2 << G4end << 105     G4cout << " ReadFile " << filename2 << G4endl;
199                                                << 106     G4int returnvalue = 0; size_t rflag = 0;
200   G4int returnvalue = 0;                       << 107     char * buffer = new char[LINEBUFFSIZE];
201   size_t rflag = 0;                            << 108 
202   char* buffer = new char[LINEBUFFSIZE];       << 109     fImplicitEndian = false;
203                                                << 110     fLittleEndian = true;
204   fImplicitEndian = false;                     << 111 
205   fLittleEndian = true;                        << 112     rflag = std::fread( buffer, 1, 128, dicom ); // The first 128 bytes
206                                                << 113                                                  //are not important
207   rflag = std::fread(buffer, 1, 128, dicom);   << 114                                                  // Reads the "DICOM" letters
208                                                << 115     rflag = std::fread( buffer, 1, 4, dicom );
209                                                << 116     // if there is no preamble, the FILE pointer is rewinded.
210   rflag = std::fread(buffer, 1, 4, dicom);     << 117     if(std::strncmp("DICM", buffer, 4) != 0) {
211   // if there is no preamble, the FILE pointer << 118         std::fseek(dicom, 0, SEEK_SET);
212   if (std::strncmp("DICM", buffer, 4) != 0) {  << 119         fImplicitEndian = true;
213     std::fseek(dicom, 0, SEEK_SET);            << 120     }
214     fImplicitEndian = true;                    << 121 
215   }                                            << 122     short readGroupId;    // identify the kind of input data
216                                                << 123     short readElementId;  // identify a particular type information
217   short readGroupId;  // identify the kind of  << 124     short elementLength2; // deal with element length in 2 bytes
218   short readElementId;  // identify a particul << 125                           //unsigned int elementLength4; 
219   short elementLength2;  // deal with element  << 126     // deal with element length in 4 bytes
220                          // unsigned int eleme << 127     unsigned long elementLength4; // deal with element length in 4 bytes
221   // deal with element length in 4 bytes       << 128 
222   unsigned long elementLength4;  // deal with  << 129     char * data = new char[DATABUFFSIZE];
223                                                << 130 
224   char* data = new char[DATABUFFSIZE];         << 131     // Read information up to the pixel data
225                                                << 132     while(true) {
226   // Read information up to the pixel data     << 133       
227   while (true) {                               << 134       //Reading groups and elements :
228     // Reading groups and elements :           << 135       readGroupId = 0;
229     readGroupId = 0;                           << 136       readElementId = 0;
230     readElementId = 0;                         << 137       // group ID
231     // group ID                                << 138       rflag = std::fread(buffer, 2, 1, dicom);
232     rflag = std::fread(buffer, 2, 1, dicom);   << 139       GetValue(buffer, readGroupId);
233     GetValue(buffer, readGroupId);             << 140       // element ID
234     // element ID                              << 141       rflag = std::fread(buffer, 2, 1, dicom);
235     rflag = std::fread(buffer, 2, 1, dicom);   << 142       GetValue(buffer, readElementId);
236     GetValue(buffer, readElementId);           << 143       
237                                                << 144       // Creating a tag to be identified afterward
238     // Creating a tag to be identified afterwa << 145       G4int tagDictionary = readGroupId*0x10000 + readElementId;
239     G4int tagDictionary = readGroupId * 0x1000 << 146       
240                                                << 147       // beginning of the pixels
241     // beginning of the pixels                 << 148       if(tagDictionary == 0x7FE00010) {
242     if (tagDictionary == 0x7FE00010) {         << 149         // Folling 2 fread's are modifications to original DICOM example 
243       // Following 2 fread's are modifications << 150         //(Jonathan Madsen)
244       // (Jonathan Madsen)                     << 151         rflag = std::fread(buffer,2,1,dicom);   // Reserved 2 bytes
245       if (!fImplicitEndian) rflag = std::fread << 152         // (not used for pixels)
246       // (not used for pixels)                 << 153         rflag = std::fread(buffer,4,1,dicom);   // Element Length  
247       rflag = std::fread(buffer, 4, 1, dicom); << 154         // (not used for pixels)
248       // (not used for pixels)                 << 155         break;      // Exit to ReadImageData()
249       break;  // Exit to ReadImageData()       << 156       }
250     }                                          << 157       
251                                                << 158       // VR or element length
252     // VR or element length                    << 159       rflag = std::fread(buffer,2,1,dicom);
253     rflag = std::fread(buffer, 2, 1, dicom);   << 160       GetValue(buffer, elementLength2);
254     GetValue(buffer, elementLength2);          << 161       
255                                                << 162       // If value representation (VR) is OB, OW, SQ, UN, added OF and UT
256     // If value representation (VR) is OB, OW, << 163       //the next length is 32 bits
257     // the next length is 32 bits              << 164       if((elementLength2 == 0x424f ||     // "OB"
258     if ((elementLength2 == 0x424f ||  // "OB"  << 165         elementLength2 == 0x574f ||     // "OW"
259          elementLength2 == 0x574f ||  // "OW"  << 166         elementLength2 == 0x464f ||     // "OF"
260          elementLength2 == 0x464f ||  // "OF"  << 167         elementLength2 == 0x5455 ||     // "UT"
261          elementLength2 == 0x5455 ||  // "UT"  << 168         elementLength2 == 0x5153 ||     // "SQ"
262          elementLength2 == 0x5153 ||  // "SQ"  << 169         elementLength2 == 0x4e55) &&    // "UN"
263          elementLength2 == 0x4e55)             << 170        !fImplicitEndian ) {             // explicit VR
264         &&  // "UN"                            << 171       
265         !fImplicitEndian)                      << 172       rflag = std::fread(buffer, 2, 1, dicom); // Skip 2 reserved bytes
266     {  // explicit VR                          << 173       
267                                                << 
268       rflag = std::fread(buffer, 2, 1, dicom); << 
269                                                << 
270       // element length                           174       // element length
271       rflag = std::fread(buffer, 4, 1, dicom);    175       rflag = std::fread(buffer, 4, 1, dicom);
272       GetValue(buffer, elementLength4);           176       GetValue(buffer, elementLength4);
273                                                << 177       
274       if (elementLength2 == 0x5153) {          << 178       if(elementLength2 == 0x5153)
275         if (elementLength4 == 0xFFFFFFFF) {    << 179         {
276           read_undefined_nested(dicom);        << 180           if(elementLength4 == 0xFFFFFFFF)
277           elementLength4 = 0;                  << 181             {
278         }                                      << 182             read_undefined_nested( dicom );
279         else {                                 << 183             elementLength4=0;
280           if (read_defined_nested(dicom, eleme << 184             }  else{
281             G4Exception("DicomHandler::ReadFil << 185             if(read_defined_nested( dicom, elementLength4 )==0){
282                         "Function read_defined << 186             G4Exception("DicomHandler::ReadFile",
                                                   >> 187                       "DICOM001",
                                                   >> 188                       FatalException,
                                                   >> 189                       "Function read_defined_nested() failed!");
                                                   >> 190             }
283           }                                       191           }
284         }                                      << 192         } else  {
285       }                                        << 
286       else {                                   << 
287         // Reading the information with data      193         // Reading the information with data
288         rflag = std::fread(data, elementLength << 194         rflag = std::fread(data, elementLength4,1,dicom);
289       }                                           195       }
290     }                                          << 196       
291     else {                                     << 197       }  else {
                                                   >> 198 
292       //  explicit with VR different than prev    199       //  explicit with VR different than previous ones
293       if (!fImplicitEndian || readGroupId == 2 << 200       if(!fImplicitEndian || readGroupId == 2) {
294         // G4cout << "Reading  DICOM files wit << 201         
295         //  element length (2 bytes)           << 202         //G4cout << "Reading  DICOM files with Explicit VR"<< G4endl;
                                                   >> 203         // element length (2 bytes)
296         rflag = std::fread(buffer, 2, 1, dicom    204         rflag = std::fread(buffer, 2, 1, dicom);
297         GetValue(buffer, elementLength2);         205         GetValue(buffer, elementLength2);
298         elementLength4 = elementLength2;          206         elementLength4 = elementLength2;
299                                                << 207         
300         rflag = std::fread(data, elementLength    208         rflag = std::fread(data, elementLength4, 1, dicom);
301       }                                        << 209         
302       else {  // Implicit VR                   << 210       } else {                                  // Implicit VR
303                                                << 211         
304         // G4cout << "Reading  DICOM files wit << 212         //G4cout << "Reading  DICOM files with Implicit VR"<< G4endl;
305                                                << 213         
306         // element length (4 bytes)               214         // element length (4 bytes)
307         if (std::fseek(dicom, -2, SEEK_CUR) != << 215         if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
308           G4Exception("DicomHandler::ReadFile" << 216           G4Exception("DicomHandler::ReadFile",
                                                   >> 217                   "DICOM001",
                                                   >> 218                   FatalException,
                                                   >> 219                   "fseek failed");
309         }                                         220         }
310                                                << 221     
311         rflag = std::fread(buffer, 4, 1, dicom    222         rflag = std::fread(buffer, 4, 1, dicom);
312         GetValue(buffer, elementLength4);         223         GetValue(buffer, elementLength4);
313                                                << 224         
314         // G4cout <<  std::hex<< elementLength << 225         //G4cout <<  std::hex<< elementLength4 << G4endl;
315                                                << 226         
316         if (elementLength4 == 0xFFFFFFFF) {    << 227         if(elementLength4 == 0xFFFFFFFF)
317           read_undefined_nested(dicom);        << 228           {
318           elementLength4 = 0;                  << 229             read_undefined_nested(dicom);
319         }                                      << 230             elementLength4=0;
320         else {                                 << 231           }  else{
321           rflag = std::fread(data, elementLeng    232           rflag = std::fread(data, elementLength4, 1, dicom);
322         }                                         233         }
                                                   >> 234         
                                                   >> 235       }
323       }                                           236       }
                                                   >> 237       
                                                   >> 238       // NULL termination
                                                   >> 239       data[elementLength4] = '\0';
                                                   >> 240       
                                                   >> 241       // analyzing information
                                                   >> 242       GetInformation(tagDictionary, data);
324     }                                             243     }
                                                   >> 244     
                                                   >> 245     G4String fnameG4DCM = G4String(filename2) + ".g4dcm";
325                                                   246 
326     // NULL termination                        << 247     // Perform functions originally written straight to file
327     data[elementLength4] = '\0';               << 248     DicomPhantomZSliceHeader* zslice = new DicomPhantomZSliceHeader(fnameG4DCM);
328                                                << 
329     // analyzing information                   << 
330     GetInformation(tagDictionary, data);       << 
331   }                                            << 
332                                                << 
333   G4String fnameG4DCM = G4String(filename2) +  << 
334                                                   249 
335   // Perform functions originally written stra << 250     std::map<G4float,G4String>::const_iterator ite;
336   DicomPhantomZSliceHeader* zslice = new Dicom << 251     for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ++ite){
337   for (auto ite = fMaterialIndices.cbegin(); i << 252       zslice->AddMaterial(ite->second);
338     zslice->AddMaterial(ite->second);          << 253     }
339   }                                            << 254     
340                                                << 255     zslice->SetNoVoxelX(fColumns/fCompression);
341   zslice->SetNoVoxelsX(fColumns / fCompression << 256     zslice->SetNoVoxelY(fRows/fCompression);
342   zslice->SetNoVoxelsY(fRows / fCompression);  << 257     zslice->SetNoVoxelZ(1);
343   zslice->SetNoVoxelsZ(1);                     << 
344                                                   258 
345   zslice->SetMinX(-fPixelSpacingX * fColumns / << 259     zslice->SetMinX(-fPixelSpacingX*fColumns/2.);
346   zslice->SetMaxX(fPixelSpacingX * fColumns /  << 260     zslice->SetMaxX(fPixelSpacingX*fColumns/2.);
347                                                   261 
348   zslice->SetMinY(-fPixelSpacingY * fRows / 2. << 262     zslice->SetMinY(-fPixelSpacingY*fRows/2.);
349   zslice->SetMaxY(fPixelSpacingY * fRows / 2.) << 263     zslice->SetMaxY(fPixelSpacingY*fRows/2.);
350                                                   264 
351   zslice->SetMinZ(fSliceLocation - fSliceThick << 265     zslice->SetMinZ(fSliceLocation-fSliceThickness/2.);
352   zslice->SetMaxZ(fSliceLocation + fSliceThick << 266     zslice->SetMaxZ(fSliceLocation+fSliceThickness/2.);
353                                                   267 
354   //===                                        << 268     //===
355                                                   269 
356   ReadData(dicom, filename2);                  << 270     ReadData( dicom, filename2 );
357                                                   271 
358   // DEPRECIATED                               << 272     // DEPRECIATED
359   // StoreData( foutG4DCM );                   << 273     //StoreData( foutG4DCM );
360   // foutG4DCM.close();                        << 274     //foutG4DCM.close();
361                                                   275 
362   StoreData(zslice);                           << 276     StoreData( zslice );
363                                                   277 
364   // Dumped 2 file after DicomPhantomZSliceMer << 278     // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency
365   // zslice->DumpToFile();                     << 279     //zslice->DumpToFile();
366                                                   280 
367   fMergedSlices->AddZSlice(zslice);            << 281     fMergedSlices->AddZSlice(zslice);
368                                                   282 
369   //                                           << 283     //
370   delete[] buffer;                             << 284     delete [] buffer;
371   delete[] data;                               << 285     delete [] data;
372                                                   286 
373   if (rflag) return returnvalue;               << 287     if (rflag) return returnvalue;
374   return returnvalue;                          << 288     return returnvalue;
375 }                                                 289 }
376                                                   290 
377 //....oooOO0OOooo........oooOO0OOooo........oo    291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378                                                   292 
379 void DicomHandler::GetInformation(G4int& tagDi << 293 void DicomHandler::GetInformation(G4int & tagDictionary, char * data)
380 {                                                 294 {
381   if (tagDictionary == 0x00280010) {  // Numbe << 295     if(tagDictionary == 0x00280010 ) { // Number of Rows
382     GetValue(data, fRows);                     << 296         GetValue(data, fRows);
383     std::printf("[0x00280010] Rows -> %i\n", f << 297         std::printf("[0x00280010] Rows -> %i\n",fRows);
384   }                                            << 298 
385   else if (tagDictionary == 0x00280011) {  //  << 299     } else if(tagDictionary == 0x00280011 ) { // Number of fColumns
386     GetValue(data, fColumns);                  << 300         GetValue(data, fColumns);
387     std::printf("[0x00280011] Columns -> %i\n" << 301         std::printf("[0x00280011] Columns -> %i\n",fColumns);
388   }                                            << 302 
389   else if (tagDictionary == 0x00280102) {  //  << 303     } else if(tagDictionary == 0x00280102 ) { // High bits  ( not used )
390     short highBits;                            << 304         short highBits;
391     GetValue(data, highBits);                  << 305         GetValue(data, highBits);
392     std::printf("[0x00280102] High bits -> %i\ << 306         std::printf("[0x00280102] High bits -> %i\n",highBits);
393   }                                            << 307 
394   else if (tagDictionary == 0x00280100) {  //  << 308     } else if(tagDictionary == 0x00280100 ) { // Bits allocated
395     GetValue(data, fBitAllocated);             << 309         GetValue(data, fBitAllocated);
396     std::printf("[0x00280100] Bits allocated - << 310         std::printf("[0x00280100] Bits allocated -> %i\n", fBitAllocated);
397   }                                            << 311 
398   else if (tagDictionary == 0x00280101) {  //  << 312     } else if(tagDictionary == 0x00280101 ) { //  Bits stored ( not used )
399     short bitStored;                           << 313         short bitStored;
400     GetValue(data, bitStored);                 << 314         GetValue(data, bitStored);
401     std::printf("[0x00280101] Bits stored -> % << 315         std::printf("[0x00280101] Bits stored -> %i\n",bitStored);
402   }                                            << 316 
403   else if (tagDictionary == 0x00280106) {  //  << 317     } else if(tagDictionary == 0x00280106 ) { //  Min. pixel value
404     GetValue(data, fMinPixelValue);            << 318         GetValue(data, fMinPixelValue);
405     std::printf("[0x00280106] Min. pixel value << 319         std::printf("[0x00280106] Min. pixel value -> %i\n", fMinPixelValue);
406   }                                            << 320 
407   else if (tagDictionary == 0x00280107) {  //  << 321     } else if(tagDictionary == 0x00280107 ) { //  Max. pixel value
408     GetValue(data, fMaxPixelValue);            << 322         GetValue(data, fMaxPixelValue);
409     std::printf("[0x00280107] Max. pixel value << 323         std::printf("[0x00280107] Max. pixel value -> %i\n", fMaxPixelValue);
410   }                                            << 324 
411   else if (tagDictionary == 0x00281053) {  //  << 325     } else if(tagDictionary == 0x00281053) { //  Rescale slope
412     fRescaleSlope = atoi(data);                << 326         fRescaleSlope = atoi(data);
413     std::printf("[0x00281053] Rescale Slope -> << 327         std::printf("[0x00281053] Rescale Slope -> %d\n", fRescaleSlope);
414   }                                            << 328 
415   else if (tagDictionary == 0x00281052) {  //  << 329     } else if(tagDictionary == 0x00281052 ) { // Rescalse intercept
416     fRescaleIntercept = atoi(data);            << 330         fRescaleIntercept = atoi(data);
417     std::printf("[0x00281052] Rescale Intercep << 331         std::printf("[0x00281052] Rescale Intercept -> %d\n", 
418   }                                            << 332                     fRescaleIntercept );
419   else if (tagDictionary == 0x00280103) {      << 333 
420     //  Pixel representation ( functions not d << 334     } else if(tagDictionary == 0x00280103 ) {
421     fPixelRepresentation = atoi(data);  // 0:  << 335         //  Pixel representation ( functions not design to read signed bits )
422     std::printf("[0x00280103] Pixel Representa << 336       fPixelRepresentation = atoi(data); // 0: unsigned  1: signed
423     if (fPixelRepresentation == 1) {           << 337         std::printf("[0x00280103] Pixel Representation -> %i\n",
424       std::printf("### PIXEL REPRESENTATION =  << 338                     fPixelRepresentation);
425       std::printf("DICOM READING SCAN FOR UNSI << 339         if(fPixelRepresentation == 1 ) {
426       std::printf("ERROR !!!!!! -> \n");       << 340             std::printf("### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
                                                   >> 341             std::printf("DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
                                                   >> 342             std::printf("ERROR !!!!!! -> \n");
                                                   >> 343         }
                                                   >> 344 
                                                   >> 345     } else if(tagDictionary == 0x00080006 ) { //  Modality
                                                   >> 346         std::printf("[0x00080006] Modality -> %s\n", data);
                                                   >> 347 
                                                   >> 348     } else if(tagDictionary == 0x00080070 ) { //  Manufacturer
                                                   >> 349         std::printf("[0x00080070] Manufacturer -> %s\n", data);
                                                   >> 350 
                                                   >> 351     } else if(tagDictionary == 0x00080080 ) { //  Institution Name
                                                   >> 352         std::printf("[0x00080080] Institution Name -> %s\n", data);
                                                   >> 353 
                                                   >> 354     } else if(tagDictionary == 0x00080081 ) { //  Institution Address
                                                   >> 355         std::printf("[0x00080081] Institution Address -> %s\n", data);
                                                   >> 356 
                                                   >> 357     } else if(tagDictionary == 0x00081040 ) { //  Institution Department Name
                                                   >> 358         std::printf("[0x00081040] Institution Department Name -> %s\n", data);
                                                   >> 359 
                                                   >> 360     } else if(tagDictionary == 0x00081090 ) { //  Manufacturer's Model Name
                                                   >> 361         std::printf("[0x00081090] Manufacturer's Model Name -> %s\n", data);
                                                   >> 362 
                                                   >> 363     } else if(tagDictionary == 0x00181000 ) { //  Device Serial Number
                                                   >> 364         std::printf("[0x00181000] Device Serial Number -> %s\n", data);
                                                   >> 365 
                                                   >> 366     } else if(tagDictionary == 0x00080008 ) { //  Image type ( not used )
                                                   >> 367         std::printf("[0x00080008] Image Types -> %s\n", data);
                                                   >> 368 
                                                   >> 369     } else if(tagDictionary == 0x00283000 ) { //Modality LUT Sequence(not used)
                                                   >> 370         std::printf("[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
                                                   >> 371 
                                                   >> 372     } else if(tagDictionary == 0x00283002 ) { // LUT Descriptor ( not used )
                                                   >> 373         std::printf("[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
                                                   >> 374 
                                                   >> 375     } else if(tagDictionary == 0x00283003 ) { // LUT Explanation ( not used )
                                                   >> 376         std::printf("[0x00283003] LUT Explanation LO 1 -> %s\n", data);
                                                   >> 377 
                                                   >> 378     } else if(tagDictionary == 0x00283004 ) { // Modality LUT ( not used )
                                                   >> 379         std::printf("[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
                                                   >> 380 
                                                   >> 381     } else if(tagDictionary == 0x00283006 ) { // LUT Data ( not used )
                                                   >> 382         std::printf("[0x00283006] LUT Data US or SS -> %s\n", data);
                                                   >> 383 
                                                   >> 384     } else if(tagDictionary == 0x00283010 ) { // VOI LUT ( not used )
                                                   >> 385         std::printf("[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
                                                   >> 386 
                                                   >> 387     } else if(tagDictionary == 0x00280120 ) { // Pixel Padding Value (not used)
                                                   >> 388       std::printf("[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
                                                   >> 389 
                                                   >> 390     } else if(tagDictionary == 0x00280030 ) { // Pixel Spacing
                                                   >> 391         G4String datas(data);
                                                   >> 392         int iss = datas.find('\\');
                                                   >> 393         fPixelSpacingX = atof( datas.substr(0,iss).c_str() );
                                                   >> 394         fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() );
                                                   >> 395 
                                                   >> 396     } else if(tagDictionary == 0x00200037 ) { // Image Orientation ( not used )
                                                   >> 397         std::printf("[0x00200037] Image Orientation (Phantom) -> %s\n", data);
                                                   >> 398 
                                                   >> 399     } else if(tagDictionary == 0x00200032 ) { // Image Position ( not used )
                                                   >> 400         std::printf("[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
                                                   >> 401 
                                                   >> 402     } else if(tagDictionary == 0x00180050 ) { // Slice Thickness
                                                   >> 403       fSliceThickness = atof(data);
                                                   >> 404       std::printf("[0x00180050] Slice Thickness (mm) -> %f\n", fSliceThickness);
                                                   >> 405 
                                                   >> 406     } else if(tagDictionary == 0x00201041 ) { // Slice Location
                                                   >> 407       fSliceLocation = atof(data);
                                                   >> 408       std::printf("[0x00201041] Slice Location -> %f\n", fSliceLocation);
                                                   >> 409 
                                                   >> 410     } else if(tagDictionary == 0x00280004 ) { // Photometric Interpretation
                                                   >> 411       // ( not used )
                                                   >> 412         std::printf("[0x00280004] Photometric Interpretation -> %s\n", data);
                                                   >> 413 
                                                   >> 414     } else if(tagDictionary == 0x00020010) { // Endian
                                                   >> 415         if(strcmp(data, "1.2.840.10008.1.2") == 0)
                                                   >> 416             fImplicitEndian = true;
                                                   >> 417         else if(strncmp(data, "1.2.840.10008.1.2.2", 19) == 0)
                                                   >> 418             fLittleEndian = false;
                                                   >> 419         //else 1.2.840..10008.1.2.1 (explicit little endian)
                                                   >> 420 
                                                   >> 421         std::printf("[0x00020010] Endian -> %s\n", data);
427     }                                             422     }
428   }                                            << 
429   else if (tagDictionary == 0x00080006) {  //  << 
430     std::printf("[0x00080006] Modality -> %s\n << 
431   }                                            << 
432   else if (tagDictionary == 0x00080070) {  //  << 
433     std::printf("[0x00080070] Manufacturer ->  << 
434   }                                            << 
435   else if (tagDictionary == 0x00080080) {  //  << 
436     std::printf("[0x00080080] Institution Name << 
437   }                                            << 
438   else if (tagDictionary == 0x00080081) {  //  << 
439     std::printf("[0x00080081] Institution Addr << 
440   }                                            << 
441   else if (tagDictionary == 0x00081040) {  //  << 
442     std::printf("[0x00081040] Institution Depa << 
443   }                                            << 
444   else if (tagDictionary == 0x00081090) {  //  << 
445     std::printf("[0x00081090] Manufacturer's M << 
446   }                                            << 
447   else if (tagDictionary == 0x00181000) {  //  << 
448     std::printf("[0x00181000] Device Serial Nu << 
449   }                                            << 
450   else if (tagDictionary == 0x00080008) {  //  << 
451     std::printf("[0x00080008] Image Types -> % << 
452   }                                            << 
453   else if (tagDictionary == 0x00283000) {  //  << 
454     std::printf("[0x00283000] Modality LUT Seq << 
455   }                                            << 
456   else if (tagDictionary == 0x00283002) {  //  << 
457     std::printf("[0x00283002] LUT Descriptor U << 
458   }                                            << 
459   else if (tagDictionary == 0x00283003) {  //  << 
460     std::printf("[0x00283003] LUT Explanation  << 
461   }                                            << 
462   else if (tagDictionary == 0x00283004) {  //  << 
463     std::printf("[0x00283004] Modality LUT Typ << 
464   }                                            << 
465   else if (tagDictionary == 0x00283006) {  //  << 
466     std::printf("[0x00283006] LUT Data US or S << 
467   }                                            << 
468   else if (tagDictionary == 0x00283010) {  //  << 
469     std::printf("[0x00283010] VOI LUT Sequence << 
470   }                                            << 
471   else if (tagDictionary == 0x00280120) {  //  << 
472     std::printf("[0x00280120] Pixel Padding Va << 
473   }                                            << 
474   else if (tagDictionary == 0x00280030) {  //  << 
475     G4String datas(data);                      << 
476     G4int iss = G4int(datas.find('\\'));       << 
477     fPixelSpacingX = atof(datas.substr(0, iss) << 
478     fPixelSpacingY = atof(datas.substr(iss + 1 << 
479   }                                            << 
480   else if (tagDictionary == 0x00200037) {  //  << 
481     std::printf("[0x00200037] Image Orientatio << 
482   }                                            << 
483   else if (tagDictionary == 0x00200032) {  //  << 
484     std::printf("[0x00200032] Image Position ( << 
485   }                                            << 
486   else if (tagDictionary == 0x00180050) {  //  << 
487     fSliceThickness = atof(data);              << 
488     std::printf("[0x00180050] Slice Thickness  << 
489   }                                            << 
490   else if (tagDictionary == 0x00201041) {  //  << 
491     fSliceLocation = atof(data);               << 
492     std::printf("[0x00201041] Slice Location - << 
493   }                                            << 
494   else if (tagDictionary == 0x00280004) {  //  << 
495     // ( not used )                            << 
496     std::printf("[0x00280004] Photometric Inte << 
497   }                                            << 
498   else if (tagDictionary == 0x00020010) {  //  << 
499     if (strcmp(data, "1.2.840.10008.1.2") == 0 << 
500       fImplicitEndian = true;                  << 
501     else if (strncmp(data, "1.2.840.10008.1.2. << 
502       fLittleEndian = false;                   << 
503     // else 1.2.840..10008.1.2.1 (explicit lit << 
504                                                   423 
505     std::printf("[0x00020010] Endian -> %s\n", << 424     // others
506   }                                            << 425     else {
                                                   >> 426         //std::printf("[0x%x] -> %s\n", tagDictionary, data);
                                                   >> 427         ;
                                                   >> 428     }
507                                                   429 
508   // others                                    << 
509   else {                                       << 
510     // std::printf("[0x%x] -> %s\n", tagDictio << 
511     ;                                          << 
512   }                                            << 
513 }                                                 430 }
514                                                   431 
515 //....oooOO0OOooo........oooOO0OOooo........oo    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
516                                                   433 
517 void DicomHandler::StoreData(DicomPhantomZSlic    434 void DicomHandler::StoreData(DicomPhantomZSliceHeader* dcmPZSH)
518 {                                                 435 {
519   G4int mean;                                  << 436     G4int mean;
520   G4double density;                            << 437     G4double density;
521   G4bool overflow = false;                     << 438     G4bool overflow = false;
522                                                << 439 
523   if (!dcmPZSH) {                              << 440     if(!dcmPZSH) { return; }
524     return;                                    << 441 
525   }                                            << 442     dcmPZSH->SetSliceLocation(fSliceLocation);
526                                                << 443 
527   dcmPZSH->SetSliceLocation(fSliceLocation);   << 444     //----- Print indices of material
528                                                << 445     if(fCompression == 1) { // no fCompression: each pixel has a density value)
529   //----- Print indices of material            << 446         for( G4int ww = 0; ww < fRows; ww++) {
530   if (fCompression == 1) {  // no fCompression << 447             dcmPZSH->AddRow();
531     for (G4int ww = 0; ww < fRows; ++ww) {     << 448             for( G4int xx = 0; xx < fColumns; xx++) {
532       dcmPZSH->AddRow();                       << 449                 mean = fTab[ww][xx];
533       for (G4int xx = 0; xx < fColumns; ++xx)  << 450                 density = Pixel2density(mean);
534         mean = fTab[ww][xx];                   << 451                 dcmPZSH->AddValue(density);
535         density = Pixel2density(mean);         << 452                 dcmPZSH->AddMateID(GetMaterialIndex(density));
536         dcmPZSH->AddValue(density);            << 453             }
537         dcmPZSH->AddMateID(GetMaterialIndex(G4 << 
538       }                                        << 
539     }                                          << 
540   }                                            << 
541   else {                                       << 
542     // density value is the average of a squar << 
543     // fCompression*fCompression pixels        << 
544     for (G4int ww = 0; ww < fRows; ww += fComp << 
545       dcmPZSH->AddRow();                       << 
546       for (G4int xx = 0; xx < fColumns; xx +=  << 
547         overflow = false;                      << 
548         mean = 0;                              << 
549         for (G4int sumx = 0; sumx < fCompressi << 
550           for (G4int sumy = 0; sumy < fCompres << 
551             if (ww + sumy >= fRows || xx + sum << 
552             mean += fTab[ww + sumy][xx + sumx] << 
553           }                                    << 
554           if (overflow) break;                 << 
555         }                                         454         }
556         mean /= fCompression * fCompression;   << 
557                                                   455 
558         if (!overflow) {                       << 456     } else {
559           density = Pixel2density(mean);       << 457         // density value is the average of a square region of
560           dcmPZSH->AddValue(density);          << 458         // fCompression*fCompression pixels
561           dcmPZSH->AddMateID(GetMaterialIndex( << 459       for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
                                                   >> 460         dcmPZSH->AddRow();
                                                   >> 461         for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
                                                   >> 462           overflow = false;
                                                   >> 463           mean = 0;
                                                   >> 464           for(int sumx = 0; sumx < fCompression; sumx++) {
                                                   >> 465             for(int sumy = 0; sumy < fCompression; sumy++) {
                                                   >> 466               if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
                                                   >> 467               mean += fTab[ww+sumy][xx+sumx];
                                                   >> 468             }
                                                   >> 469             if(overflow) break;
                                                   >> 470           }
                                                   >> 471           mean /= fCompression*fCompression;
                                                   >> 472           
                                                   >> 473           if(!overflow) {
                                                   >> 474             density = Pixel2density(mean);
                                                   >> 475             dcmPZSH->AddValue(density);
                                                   >> 476             dcmPZSH->AddMateID(GetMaterialIndex(density));
                                                   >> 477           }
562         }                                         478         }
563       }                                           479       }
564     }                                             480     }
565   }                                            << 481     
566                                                << 482     dcmPZSH->FlipData();
567   dcmPZSH->FlipData();                         << 
568 }                                                 483 }
569                                                   484 
570 //....oooOO0OOooo........oooOO0OOooo........oo    485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
571 // This function is depreciated as it is handl << 486 // This function is depreciated as it is handled by 
572 // DicomPhantomZSliceHeader::DumpToFile           487 // DicomPhantomZSliceHeader::DumpToFile
573 void DicomHandler::StoreData(std::ofstream& fo    488 void DicomHandler::StoreData(std::ofstream& foutG4DCM)
574 {                                                 489 {
575   G4int mean;                                     490   G4int mean;
576   G4double density;                               491   G4double density;
577   G4bool overflow = false;                        492   G4bool overflow = false;
578                                                << 493   
579   //----- Print indices of material               494   //----- Print indices of material
580   if (fCompression == 1) {  // no fCompression << 495   if(fCompression == 1) { // no fCompression: each pixel has a density value)
581     for (G4int ww = 0; ww < fRows; ++ww) {     << 496     for( G4int ww = 0; ww < fRows; ww++) {
582       for (G4int xx = 0; xx < fColumns; ++xx)  << 497       for( G4int xx = 0; xx < fColumns; xx++) {
583         mean = fTab[ww][xx];                      498         mean = fTab[ww][xx];
584         density = Pixel2density(mean);            499         density = Pixel2density(mean);
585         foutG4DCM << GetMaterialIndex(G4float( << 500         foutG4DCM << GetMaterialIndex( density ) << " ";
586       }                                           501       }
587       foutG4DCM << G4endl;                        502       foutG4DCM << G4endl;
588     }                                             503     }
589   }                                            << 504     
590   else {                                       << 505   } else {
591     // density value is the average of a squar    506     // density value is the average of a square region of
592     // fCompression*fCompression pixels           507     // fCompression*fCompression pixels
593     for (G4int ww = 0; ww < fRows; ww += fComp << 508     for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
594       for (G4int xx = 0; xx < fColumns; xx +=  << 509       for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
595         overflow = false;                         510         overflow = false;
596         mean = 0;                                 511         mean = 0;
597         for (G4int sumx = 0; sumx < fCompressi << 512         for(int sumx = 0; sumx < fCompression; sumx++) {
598           for (G4int sumy = 0; sumy < fCompres << 513           for(int sumy = 0; sumy < fCompression; sumy++) {
599             if (ww + sumy >= fRows || xx + sum << 514             if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
600             mean += fTab[ww + sumy][xx + sumx] << 515             mean += fTab[ww+sumy][xx+sumx];
601           }                                       516           }
602           if (overflow) break;                 << 517           if(overflow) break;
603         }                                         518         }
604         mean /= fCompression * fCompression;   << 519         mean /= fCompression*fCompression;
605                                                << 520         
606         if (!overflow) {                       << 521         if(!overflow) {
607           density = Pixel2density(mean);          522           density = Pixel2density(mean);
608           foutG4DCM << GetMaterialIndex(G4floa << 523           foutG4DCM << GetMaterialIndex( density ) << " ";
609         }                                         524         }
610       }                                           525       }
611       foutG4DCM << G4endl;                        526       foutG4DCM << G4endl;
612     }                                             527     }
                                                   >> 528     
613   }                                               529   }
614                                                << 530   
615   //----- Print densities                         531   //----- Print densities
616   if (fCompression == 1) {  // no fCompression << 532   if(fCompression == 1) { // no fCompression: each pixel has a density value)
617     for (G4int ww = 0; ww < fRows; ww++) {     << 533     for( G4int ww = 0; ww < fRows; ww++) {
618       for (G4int xx = 0; xx < fColumns; xx++)  << 534       for( G4int xx = 0; xx < fColumns; xx++) {
619         mean = fTab[ww][xx];                      535         mean = fTab[ww][xx];
620         density = Pixel2density(mean);            536         density = Pixel2density(mean);
621         foutG4DCM << density << " ";              537         foutG4DCM << density << " ";
622         if (xx % 8 == 3) foutG4DCM << G4endl;  << 538         if( xx%8 == 3 ) foutG4DCM << G4endl; // just for nicer reading
623       }                                           539       }
624     }                                             540     }
625   }                                            << 541     
626   else {                                       << 542   } else {
627     // density value is the average of a squar    543     // density value is the average of a square region of
628     // fCompression*fCompression pixels           544     // fCompression*fCompression pixels
629     for (G4int ww = 0; ww < fRows; ww += fComp << 545     for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
630       for (G4int xx = 0; xx < fColumns; xx +=  << 546       for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
631         overflow = false;                         547         overflow = false;
632         mean = 0;                                 548         mean = 0;
633         for (G4int sumx = 0; sumx < fCompressi << 549         for(int sumx = 0; sumx < fCompression; sumx++) {
634           for (G4int sumy = 0; sumy < fCompres << 550           for(int sumy = 0; sumy < fCompression; sumy++) {
635             if (ww + sumy >= fRows || xx + sum << 551             if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
636             mean += fTab[ww + sumy][xx + sumx] << 552             mean += fTab[ww+sumy][xx+sumx];
637           }                                       553           }
638           if (overflow) break;                 << 554           if(overflow) break;
639         }                                         555         }
640         mean /= fCompression * fCompression;   << 556         mean /= fCompression*fCompression;
641                                                << 557         
642         if (!overflow) {                       << 558         if(!overflow) {
643           density = Pixel2density(mean);          559           density = Pixel2density(mean);
644           foutG4DCM << density << " ";         << 560           foutG4DCM << density  << " ";
645           if (xx / fCompression % 8 == 3) fout << 561           if( xx/fCompression%8 == 3 ) foutG4DCM << G4endl; // just for nicer
646           // reading                              562           // reading
647         }                                         563         }
648       }                                           564       }
649     }                                             565     }
                                                   >> 566     
650   }                                               567   }
                                                   >> 568   
651 }                                                 569 }
652                                                   570 
653 //....oooOO0OOooo........oooOO0OOooo........oo << 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
654 void DicomHandler::ReadMaterialIndices(std::if << 572 void DicomHandler::ReadMaterialIndices( std::ifstream& finData)
655 {                                                 573 {
656   unsigned int nMate;                             574   unsigned int nMate;
657   G4String mateName;                              575   G4String mateName;
658   G4float densityMax;                             576   G4float densityMax;
659   finData >> nMate;                               577   finData >> nMate;
660   if (finData.eof()) return;                   << 578   if( finData.eof() ) return;
661                                                << 579   
662   G4cout << " ReadMaterialIndices " << nMate <    580   G4cout << " ReadMaterialIndices " << nMate << G4endl;
663   for (unsigned int ii = 0; ii < nMate; ++ii)  << 581   for( unsigned int ii = 0; ii < nMate; ii++ ){
664     finData >> mateName >> densityMax;            582     finData >> mateName >> densityMax;
665     fMaterialIndices[densityMax] = mateName;      583     fMaterialIndices[densityMax] = mateName;
666     //    G4cout << ii << " ReadMaterialIndice << 584     G4cout << ii << " ReadMaterialIndices " << mateName << " " 
667     //     << densityMax << G4endl;            << 585            << densityMax << G4endl;
668   }                                               586   }
                                                   >> 587   
669 }                                                 588 }
670                                                   589 
671 //....oooOO0OOooo........oooOO0OOooo........oo    590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
672                                                   591 
673 unsigned int DicomHandler::GetMaterialIndex(G4 << 592 unsigned int DicomHandler::GetMaterialIndex( G4float density )
674 {                                                 593 {
675   std::map<G4float, G4String>::const_reverse_i << 594   std::map<G4float,G4String>::reverse_iterator ite;
676   std::size_t ii = fMaterialIndices.size();    << 595   G4int ii = fMaterialIndices.size();
677                                                << 596   for( ite = fMaterialIndices.rbegin(); ite != fMaterialIndices.rend();
678   for (ite = fMaterialIndices.crbegin(); ite ! << 597        ite++, ii-- ) {
679     if (density >= (*ite).first) {             << 598     if( density >= (*ite).first ) {
680       break;                                      599       break;
681     }                                             600     }
682   }                                               601   }
683                                                << 602   //-  G4cout << " GetMaterialIndex " << density << " = " << ii << G4endl;
684   if (ii == fMaterialIndices.size()) {         << 603   
685     ii = fMaterialIndices.size() - 1;          << 604   if(static_cast<unsigned int>(ii) == fMaterialIndices.size())
686   }                                            << 605     { ii = fMaterialIndices.size()-1; }
687                                                << 606   
688   return unsigned(ii);                         << 607   return  ii;
                                                   >> 608   
689 }                                                 609 }
690                                                   610 
691 //....oooOO0OOooo........oooOO0OOooo........oo    611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
692                                                   612 
693 G4int DicomHandler::ReadData(FILE* dicom, char << 613 G4int DicomHandler::ReadData(FILE *dicom,char * filename2)
694 {                                                 614 {
695   G4int returnvalue = 0;                       << 615   G4int returnvalue = 0; size_t rflag = 0;
696   size_t rflag = 0;                            << 616   
697                                                << 617   //  READING THE PIXELS :
698   //  READING THE PIXELS                       << 618   G4int w = 0;
699                                                << 619   
700   fTab = new G4int*[fRows];                       620   fTab = new G4int*[fRows];
701   for (G4int i = 0; i < fRows; ++i) {          << 621   for ( G4int i = 0; i < fRows; i ++ ) {
702     fTab[i] = new G4int[fColumns];             << 622         fTab[i] = new G4int[fColumns];
703   }                                            << 623     }
704                                                << 624   
705   if (fBitAllocated == 8) {  // Case 8 bits :  << 625   if(fBitAllocated == 8) { // Case 8 bits :
706                                                << 626     
707     std::printf("@@@ Error! Picture != 16 bits    627     std::printf("@@@ Error! Picture != 16 bits...\n");
708     std::printf("@@@ Error! Picture != 16 bits    628     std::printf("@@@ Error! Picture != 16 bits...\n");
709     std::printf("@@@ Error! Picture != 16 bits    629     std::printf("@@@ Error! Picture != 16 bits...\n");
710                                                << 630     
711     unsigned char ch = 0;                         631     unsigned char ch = 0;
712                                                << 632     
713     for (G4int j = 0; j < fRows; ++j) {        << 633     for(G4int j = 0; j < fRows; j++) {
714       for (G4int i = 0; i < fColumns; ++i) {   << 634       for(G4int i = 0; i < fColumns; i++) {
715         rflag = std::fread(&ch, 1, 1, dicom);  << 635         w++;
716         fTab[j][i] = ch * fRescaleSlope + fRes << 636         rflag = std::fread( &ch, 1, 1, dicom);
                                                   >> 637         fTab[j][i] = ch*fRescaleSlope + fRescaleIntercept;
717       }                                           638       }
718     }                                             639     }
719     returnvalue = 1;                              640     returnvalue = 1;
720   }                                            << 641     
721   else {  //  from 12 to 16 bits :             << 642   } else { //  from 12 to 16 bits :
722     char sbuff[2];                                643     char sbuff[2];
723     short pixel;                                  644     short pixel;
724     for (G4int j = 0; j < fRows; ++j) {        << 645     for( G4int j = 0; j < fRows; j++) {
725       for (G4int i = 0; i < fColumns; ++i) {   << 646       for( G4int i = 0; i < fColumns; i++) {
                                                   >> 647         w++;
726         rflag = std::fread(sbuff, 2, 1, dicom)    648         rflag = std::fread(sbuff, 2, 1, dicom);
727         GetValue(sbuff, pixel);                   649         GetValue(sbuff, pixel);
728         fTab[j][i] = pixel * fRescaleSlope + f << 650         fTab[j][i] = pixel*fRescaleSlope + fRescaleIntercept;
729       }                                           651       }
730     }                                             652     }
731   }                                               653   }
732                                                << 654   
733   // Creation of .g4 files wich contains avera    655   // Creation of .g4 files wich contains averaged density data
734   G4String nameProcessed = filename2 + G4Strin << 656   char * nameProcessed = new char[FILENAMESIZE];
735   FILE* fileOut = std::fopen(nameProcessed.c_s << 657   FILE* fileOut;
736                                                << 658   
737   G4cout << "### Writing of " << nameProcessed << 659   std::sprintf(nameProcessed,"%s.g4dcmb",filename2);
738                                                << 660   fileOut = std::fopen(nameProcessed,"w+b");
                                                   >> 661   std::printf("### Writing of %s ###\n",nameProcessed);
                                                   >> 662   
739   unsigned int nMate = fMaterialIndices.size()    663   unsigned int nMate = fMaterialIndices.size();
740   rflag = std::fwrite(&nMate, sizeof(unsigned     664   rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
741   //--- Write materials                           665   //--- Write materials
742   for (auto ite = fMaterialIndices.cbegin(); i << 666   std::map<G4float,G4String>::const_iterator ite;
                                                   >> 667   for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++ ){
743     G4String mateName = (*ite).second;            668     G4String mateName = (*ite).second;
744     for (std::size_t ii = (*ite).second.length << 669     for( G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
745       mateName += " ";                            670       mateName += " ";
746     }  // mateName = const_cast<char*>(((*ite) << 671     }         //mateName = const_cast<char*>(((*ite).second).c_str());
747                                                << 672     
748     const char* mateNameC = mateName.c_str();     673     const char* mateNameC = mateName.c_str();
749     rflag = std::fwrite(mateNameC, sizeof(char << 674     rflag = std::fwrite(mateNameC, sizeof(char),40, fileOut);
750   }                                               675   }
751                                                << 676   
752   unsigned int fRowsC = fRows / fCompression;  << 677   unsigned int fRowsC = fRows/fCompression;
753   unsigned int fColumnsC = fColumns / fCompres << 678   unsigned int fColumnsC = fColumns/fCompression;
754   unsigned int planesC = 1;                       679   unsigned int planesC = 1;
755   G4float pixelLocationXM = -G4float(fPixelSpa << 680   G4float pixelLocationXM = -fPixelSpacingX*fColumns/2.;
756   G4float pixelLocationXP = G4float(fPixelSpac << 681   G4float pixelLocationXP = fPixelSpacingX*fColumns/2.;
757   G4float pixelLocationYM = -G4float(fPixelSpa << 682   G4float pixelLocationYM = -fPixelSpacingY*fRows/2.;
758   G4float pixelLocationYP = G4float(fPixelSpac << 683   G4float pixelLocationYP = fPixelSpacingY*fRows/2.;
759   G4float fSliceLocationZM = G4float(fSliceLoc << 684   G4float fSliceLocationZM = fSliceLocation-fSliceThickness/2.;
760   G4float fSliceLocationZP = G4float(fSliceLoc << 685   G4float fSliceLocationZP = fSliceLocation+fSliceThickness/2.;
761   //--- Write number of voxels (assume only on    686   //--- Write number of voxels (assume only one voxel in Z)
762   rflag = std::fwrite(&fRowsC, sizeof(unsigned    687   rflag = std::fwrite(&fRowsC, sizeof(unsigned int), 1, fileOut);
763   rflag = std::fwrite(&fColumnsC, sizeof(unsig    688   rflag = std::fwrite(&fColumnsC, sizeof(unsigned int), 1, fileOut);
764   rflag = std::fwrite(&planesC, sizeof(unsigne    689   rflag = std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut);
765   //--- Write minimum and maximum extensions      690   //--- Write minimum and maximum extensions
766   rflag = std::fwrite(&pixelLocationXM, sizeof    691   rflag = std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut);
767   rflag = std::fwrite(&pixelLocationXP, sizeof    692   rflag = std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut);
768   rflag = std::fwrite(&pixelLocationYM, sizeof    693   rflag = std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut);
769   rflag = std::fwrite(&pixelLocationYP, sizeof    694   rflag = std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut);
770   rflag = std::fwrite(&fSliceLocationZM, sizeo    695   rflag = std::fwrite(&fSliceLocationZM, sizeof(G4float), 1, fileOut);
771   rflag = std::fwrite(&fSliceLocationZP, sizeo    696   rflag = std::fwrite(&fSliceLocationZP, sizeof(G4float), 1, fileOut);
772   // rflag = std::fwrite(&fCompression, sizeof    697   // rflag = std::fwrite(&fCompression, sizeof(unsigned int), 1, fileOut);
773                                                << 698   
774   std::printf("%8i   %8i\n", fRows, fColumns); << 699   std::printf("%8i   %8i\n",fRows,fColumns);
775   std::printf("%8f   %8f\n", fPixelSpacingX, f << 700   std::printf("%8f   %8f\n",fPixelSpacingX,fPixelSpacingY);
776   std::printf("%8f\n", fSliceThickness);          701   std::printf("%8f\n", fSliceThickness);
777   std::printf("%8f\n", fSliceLocation);           702   std::printf("%8f\n", fSliceLocation);
778   std::printf("%8i\n", fCompression);             703   std::printf("%8i\n", fCompression);
779                                                << 704   
780   G4int compSize = fCompression;                  705   G4int compSize = fCompression;
781   G4int mean;                                     706   G4int mean;
782   G4float density;                                707   G4float density;
783   G4bool overflow = false;                        708   G4bool overflow = false;
784                                                << 709   
785   //----- Write index of material for each pix    710   //----- Write index of material for each pixel
786   if (compSize == 1) {  // no fCompression: ea << 711   if(compSize == 1) { // no fCompression: each pixel has a density value)
787     for (G4int ww = 0; ww < fRows; ++ww) {     << 712     for( G4int ww = 0; ww < fRows; ww++) {
788       for (G4int xx = 0; xx < fColumns; ++xx)  << 713       for( G4int xx = 0; xx < fColumns; xx++) {
789         mean = fTab[ww][xx];                      714         mean = fTab[ww][xx];
790         density = Pixel2density(mean);            715         density = Pixel2density(mean);
791         unsigned int mateID = GetMaterialIndex << 716         unsigned int mateID = GetMaterialIndex( density );
792         rflag = std::fwrite(&mateID, sizeof(un    717         rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
793       }                                           718       }
794     }                                             719     }
795   }                                            << 720     
796   else {                                       << 721   } else {
797     // density value is the average of a squar    722     // density value is the average of a square region of
798     // fCompression*fCompression pixels           723     // fCompression*fCompression pixels
799     for (G4int ww = 0; ww < fRows; ww += compS << 724     for(G4int ww = 0; ww < fRows ;ww += compSize ) {
800       for (G4int xx = 0; xx < fColumns; xx +=  << 725       for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
801         overflow = false;                         726         overflow = false;
802         mean = 0;                                 727         mean = 0;
803         for (G4int sumx = 0; sumx < compSize;  << 728         for(int sumx = 0; sumx < compSize; sumx++) {
804           for (G4int sumy = 0; sumy < compSize << 729           for(int sumy = 0; sumy < compSize; sumy++) {
805             if (ww + sumy >= fRows || xx + sum << 730             if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
806             mean += fTab[ww + sumy][xx + sumx] << 731             mean += fTab[ww+sumy][xx+sumx];
807           }                                       732           }
808           if (overflow) break;                 << 733           if(overflow) break;
809         }                                         734         }
810         mean /= compSize * compSize;           << 735         mean /= compSize*compSize;
811                                                << 736         
812         if (!overflow) {                       << 737         if(!overflow) {
813           density = Pixel2density(mean);          738           density = Pixel2density(mean);
814           unsigned int mateID = GetMaterialInd << 739           unsigned int mateID = GetMaterialIndex( density );
815           rflag = std::fwrite(&mateID, sizeof(    740           rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
816         }                                         741         }
817       }                                           742       }
                                                   >> 743       
818     }                                             744     }
819   }                                               745   }
820                                                << 746   
821   //----- Write density for each pixel            747   //----- Write density for each pixel
822   if (compSize == 1) {  // no fCompression: ea << 748   if(compSize == 1) { // no fCompression: each pixel has a density value)
823     for (G4int ww = 0; ww < fRows; ++ww) {     << 749     for( G4int ww = 0; ww < fRows; ww++) {
824       for (G4int xx = 0; xx < fColumns; ++xx)  << 750       for( G4int xx = 0; xx < fColumns; xx++) {
825         mean = fTab[ww][xx];                      751         mean = fTab[ww][xx];
826         density = Pixel2density(mean);            752         density = Pixel2density(mean);
827         rflag = std::fwrite(&density, sizeof(G    753         rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
828       }                                           754       }
829     }                                             755     }
830   }                                            << 756     
831   else {                                       << 757   } else {
832     // density value is the average of a squar    758     // density value is the average of a square region of
833     // fCompression*fCompression pixels           759     // fCompression*fCompression pixels
834     for (G4int ww = 0; ww < fRows; ww += compS << 760     for(G4int ww = 0; ww < fRows ;ww += compSize ) {
835       for (G4int xx = 0; xx < fColumns; xx +=  << 761       for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
836         overflow = false;                         762         overflow = false;
837         mean = 0;                                 763         mean = 0;
838         for (G4int sumx = 0; sumx < compSize;  << 764         for(int sumx = 0; sumx < compSize; sumx++) {
839           for (G4int sumy = 0; sumy < compSize << 765           for(int sumy = 0; sumy < compSize; sumy++) {
840             if (ww + sumy >= fRows || xx + sum << 766             if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
841             mean += fTab[ww + sumy][xx + sumx] << 767             mean += fTab[ww+sumy][xx+sumx];
842           }                                       768           }
843           if (overflow) break;                 << 769           if(overflow) break;
844         }                                         770         }
845         mean /= compSize * compSize;           << 771         mean /= compSize*compSize;
846                                                << 772         
847         if (!overflow) {                       << 773         if(!overflow) {
848           density = Pixel2density(mean);          774           density = Pixel2density(mean);
849           rflag = std::fwrite(&density, sizeof    775           rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
850         }                                         776         }
851       }                                           777       }
                                                   >> 778       
852     }                                             779     }
853   }                                               780   }
854                                                << 781   
855   rflag = std::fclose(fileOut);                   782   rflag = std::fclose(fileOut);
856                                                << 783   
857   delete[] nameProcessed;                      << 784   delete [] nameProcessed;
858                                                << 785   
859   /*    for ( G4int i = 0; i < fRows; i ++ ) {    786   /*    for ( G4int i = 0; i < fRows; i ++ ) {
860         delete [] fTab[i];                        787         delete [] fTab[i];
861         }                                         788         }
862      delete [] fTab;                              789      delete [] fTab;
863   */                                              790   */
864                                                << 791   
865   if (rflag) return returnvalue;                  792   if (rflag) return returnvalue;
866   return returnvalue;                             793   return returnvalue;
867 }                                                 794 }
868                                                   795 
869 //....oooOO0OOooo........oooOO0OOooo........oo << 796 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
870                                                << 
871 // DICOM HEAD does not use the calibration cur << 
872                                                << 
873 #ifdef DICOM_USE_HEAD                          << 
874 void DicomHandler::ReadCalibration()           << 
875 {                                              << 
876   fNbrequali = 0;                              << 
877   fReadCalibration = false;                    << 
878   G4cout << "No calibration curve for the DICO << 
879 }                                              << 
880 #else                                          << 
881 // Separated out of Pixel2density                 797 // Separated out of Pixel2density
882 // No need to read in same calibration EVERY t    798 // No need to read in same calibration EVERY time
883 // Increases the speed of reading file by seve    799 // Increases the speed of reading file by several orders of magnitude
884                                                << 
885 void DicomHandler::ReadCalibration()              800 void DicomHandler::ReadCalibration()
886 {                                                 801 {
887   fNbrequali = 0;                                 802   fNbrequali = 0;
                                                   >> 803   
888   // CT2Density.dat contains the calibration c    804   // CT2Density.dat contains the calibration curve to convert CT (Hounsfield)
889   // number to physical density                   805   // number to physical density
890   std::ifstream calibration(fCt2DensityFile.c_    806   std::ifstream calibration(fCt2DensityFile.c_str());
891   calibration >> fNbrequali;                      807   calibration >> fNbrequali;
                                                   >> 808   
892   fValueDensity = new G4double[fNbrequali];       809   fValueDensity = new G4double[fNbrequali];
893   fValueCT = new G4double[fNbrequali];            810   fValueCT = new G4double[fNbrequali];
894                                                << 811   
895   if (!calibration) {                          << 812   if(!calibration) {
896     G4Exception("DicomHandler::ReadFile", "DIC << 813     G4Exception("DicomHandler::ReadFile",
                                                   >> 814                 "DICOM001",
                                                   >> 815                 FatalException,
897                 "@@@ No value to transform pix    816                 "@@@ No value to transform pixels in density!");
898   }                                            << 817     
899   else {  // calibration was successfully open << 818   } else { // calibration was successfully opened
900     for (G4int i = 0; i < fNbrequali; ++i) {   << 819     for(G4int i = 0; i < fNbrequali; i++) { // Loop to store all the 
                                                   >> 820       //pts in CT2Density.dat
901       calibration >> fValueCT[i] >> fValueDens    821       calibration >> fValueCT[i] >> fValueDensity[i];
902     }                                             822     }
903   }                                               823   }
904   calibration.close();                            824   calibration.close();
                                                   >> 825   
905   fReadCalibration = true;                        826   fReadCalibration = true;
906 }                                                 827 }
907 #endif                                         << 
908                                                   828 
909 #ifdef DICOM_USE_HEAD                          << 829 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
910 G4float DicomHandler::Pixel2density(G4int pixe << 
911 {                                              << 
912   G4float density = -1;                        << 
913                                                << 
914   // Air                                       << 
915   if (pixel == -998.) density = 0.001290;      << 
916   // Soft Tissue                               << 
917   else if (pixel == 24.)                       << 
918     density = 1.00;                            << 
919   // Brain                                     << 
920   else if (pixel == 52.)                       << 
921     density = 1.03;                            << 
922   // Spinal disc                               << 
923   else if (pixel == 92.)                       << 
924     density = 1.10;                            << 
925   // Trabecular bone                           << 
926   else if (pixel == 197.)                      << 
927     density = 1.18;                            << 
928   // Cortical Bone                             << 
929   else if (pixel == 923.)                      << 
930     density = 1.85;                            << 
931   // Tooth dentine                             << 
932   else if (pixel == 1280.)                     << 
933     density = 2.14;                            << 
934   // Tooth enamel                              << 
935   else if (pixel == 2310.)                     << 
936     density = 2.89;                            << 
937                                                << 
938   return density;                              << 
939 }                                              << 
940                                                << 
941 #else                                          << 
942 //....oooOO0OOooo........oooOO0OOooo........oo << 
943                                                   830 
944 G4float DicomHandler::Pixel2density(G4int pixe    831 G4float DicomHandler::Pixel2density(G4int pixel)
945 {                                                 832 {
946   if (!fReadCalibration) {                     << 833   if(!fReadCalibration) { ReadCalibration(); }
947     ReadCalibration();                         << 834   
948   }                                            << 
949                                                << 
950   G4float density = -1.;                          835   G4float density = -1.;
951   G4double deltaCT = 0;                           836   G4double deltaCT = 0;
952   G4double deltaDensity = 0;                      837   G4double deltaDensity = 0;
953                                                << 838   
954   for (G4int j = 1; j < fNbrequali; ++j) {     << 839   
955     if (pixel >= fValueCT[j - 1] && pixel < fV << 840   for(G4int j = 1; j < fNbrequali; j++) {
956       deltaCT = fValueCT[j] - fValueCT[j - 1]; << 841     if( pixel >= fValueCT[j-1] && pixel < fValueCT[j]) {
957       deltaDensity = fValueDensity[j] - fValue << 842       
958                                                << 843       deltaCT = fValueCT[j] - fValueCT[j-1];
                                                   >> 844       deltaDensity = fValueDensity[j] - fValueDensity[j-1];
                                                   >> 845       
959       // interpolating linearly                   846       // interpolating linearly
960       density = G4float(fValueDensity[j] - ((f << 847       density = fValueDensity[j] - ((fValueCT[j] - pixel)*deltaDensity/deltaCT );
961       break;                                      848       break;
962     }                                             849     }
963   }                                               850   }
964                                                << 851   
965   if (density < 0.) {                          << 852   if(density < 0.) {
966     std::printf(                               << 853     std::printf("@@@ Error density = %f && Pixel = %i \
967       "@@@ Error density = %f && Pixel = %i \  << 854       (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
968       (0x%x) && deltaDensity/deltaCT = %f\n",  << 855                 deltaDensity/deltaCT);
969       density, pixel, pixel, deltaDensity / de << 
970   }                                               856   }
971                                                << 857   
972   return density;                                 858   return density;
973 }                                                 859 }
974 #endif                                         << 860 
975 //....oooOO0OOooo........oooOO0OOooo........oo    861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
976                                                   862 
977 void DicomHandler::CheckFileFormat()              863 void DicomHandler::CheckFileFormat()
978 {                                                 864 {
979   std::ifstream checkData(fDriverFile.c_str())    865   std::ifstream checkData(fDriverFile.c_str());
980   char* oneLine = new char[128];               << 866   char * oneLine = new char[128];
981                                                << 867   
982   if (!(checkData.is_open())) {  // Check exis << 868   if(!(checkData.is_open())) { //Check existance of Data.dat
983                                                << 869     
984     G4String message = "\nDicomG4 needs Data.d << 870     G4String message = 
                                                   >> 871       "\nDicomG4 needs Data.dat (or another driver file specified";
985     message += " in command line):\n";            872     message += " in command line):\n";
986     message += "\tFirst line: number of image     873     message += "\tFirst line: number of image pixel for a voxel (G4Box)\n";
987     message += "\tSecond line: number of image    874     message += "\tSecond line: number of images (CT slices) to read\n";
988     message += "\tEach following line contains    875     message += "\tEach following line contains the name of a Dicom image";
989     message += " except for the .dcm extension    876     message += " except for the .dcm extension";
990     G4Exception("DicomHandler::ReadFile", "DIC << 877     G4Exception("DicomHandler::ReadFile",
                                                   >> 878                 "DICOM001",
                                                   >> 879                 FatalException,
                                                   >> 880                 message.c_str());
991   }                                               881   }
992                                                << 882   
993   checkData >> fCompression;                      883   checkData >> fCompression;
994   checkData >> fNFiles;                           884   checkData >> fNFiles;
995   G4String oneName;                               885   G4String oneName;
996   checkData.getline(oneLine, 100);             << 886   checkData.getline(oneLine,100);
997   std::ifstream testExistence;                    887   std::ifstream testExistence;
998   G4bool existAlready = true;                     888   G4bool existAlready = true;
999   for (G4int rep = 0; rep < fNFiles; ++rep) {  << 889   for(G4int rep = 0; rep < fNFiles; rep++) {
1000     checkData.getline(oneLine, 100);          << 890     checkData.getline(oneLine,100);
1001     oneName = oneLine;                           891     oneName = oneLine;
1002     oneName += ".g4dcm";  // create dicomFile << 892     oneName += ".g4dcm"; // create dicomFile.g4dcm
1003     G4cout << fNFiles << " test file " << one    893     G4cout << fNFiles << " test file " << oneName << G4endl;
1004     testExistence.open(oneName.data());          894     testExistence.open(oneName.data());
1005     if (!(testExistence.is_open())) {         << 895     if(!(testExistence.is_open())) {
1006       existAlready = false;                      896       existAlready = false;
1007       testExistence.clear();                     897       testExistence.clear();
1008       testExistence.close();                     898       testExistence.close();
1009     }                                            899     }
1010     testExistence.clear();                       900     testExistence.clear();
1011     testExistence.close();                       901     testExistence.close();
1012   }                                              902   }
1013                                               << 903   
1014   ReadMaterialIndices(checkData);             << 904   ReadMaterialIndices( checkData );
1015                                               << 905   
1016   checkData.close();                             906   checkData.close();
1017   delete[] oneLine;                           << 907   delete [] oneLine;
1018                                               << 908   
1019   if (existAlready == false) {  // The files  << 909   if( existAlready == false  ) { // The files *.g4dcm have to be created
1020                                               << 910     
1021     G4cout << "\nAll the necessary images wer    911     G4cout << "\nAll the necessary images were not found in processed form "
1022            << ", starting with .dcm images\n"    912            << ", starting with .dcm images\n";
1023                                               << 913     
1024     FILE* dicom;                              << 914     FILE * dicom;
1025     FILE* lecturePref;                        << 915     FILE * lecturePref;
1026     char* fCompressionc = new char[LINEBUFFSI << 916     char * fCompressionc = new char[LINEBUFFSIZE];
1027     char* maxc = new char[LINEBUFFSIZE];      << 917     char * maxc = new char[LINEBUFFSIZE];
1028     // char name[300], inputFile[300];        << 918     //char name[300], inputFile[300];
1029     char* inputFile = new char[FILENAMESIZE]; << 919     char * name = new char[FILENAMESIZE];
                                                   >> 920     char * inputFile = new char[FILENAMESIZE];
1030     G4int rflag;                                 921     G4int rflag;
1031     lecturePref = std::fopen(fDriverFile.c_st << 922     
1032                                               << 923     lecturePref = std::fopen(fDriverFile.c_str(),"r");
1033     rflag = std::fscanf(lecturePref, "%s", fC << 924     rflag = std::fscanf(lecturePref,"%s",fCompressionc);
1034     fCompression = atoi(fCompressionc);          925     fCompression = atoi(fCompressionc);
1035     rflag = std::fscanf(lecturePref, "%s", ma << 926     rflag = std::fscanf(lecturePref,"%s",maxc);
1036     fNFiles = atoi(maxc);                        927     fNFiles = atoi(maxc);
1037     G4cout << " fNFiles " << fNFiles << G4end    928     G4cout << " fNFiles " << fNFiles << G4endl;
1038                                               << 929     
1039     /////////////////////                     << 930     for( G4int i = 1; i <= fNFiles; i++ ) { // Begin loop on filenames
1040                                               << 931       
1041 #ifdef DICOM_USE_HEAD                         << 932       rflag = std::fscanf(lecturePref,"%s",inputFile);
1042     for (G4int i = 1; i <= fNFiles; ++i) {  / << 933       std::sprintf(name,"%s.dcm",inputFile);
1043       rflag = std::fscanf(lecturePref, "%s",  << 934       G4cout << "check 1: " << name << G4endl;
1044       G4String path = GetDicomDataPath() + "/ << 
1045                                               << 
1046       G4String name = inputFile + G4String(". << 
1047       // Writes the results to a character st << 
1048                                               << 
1049       G4String name2 = path + name;           << 
1050       //  Open input file and give it to gest    935       //  Open input file and give it to gestion_dicom :
1051       G4cout << "### Opening " << name2 << "  << 936       std::printf("### Opening %s and reading :\n",name);
1052       dicom = std::fopen(name2.c_str(), "rb") << 937       dicom = std::fopen(name,"rb");
1053       // Reading the .dcm in two steps:          938       // Reading the .dcm in two steps:
1054       //      1.  reading the header             939       //      1.  reading the header
1055       //      2. reading the pixel data and s    940       //      2. reading the pixel data and store the density in Moyenne.dat
1056       if (dicom != 0) {                       << 941       if( dicom != 0 ) {
1057         ReadFile(dicom, inputFile);           << 942         ReadFile(dicom,inputFile);
1058       }                                       << 943       } else {
1059       else {                                  << 
1060         G4cout << "\nError opening file : " < << 
1061       }                                       << 
1062       rflag = std::fclose(dicom);             << 
1063     }                                         << 
1064 #else                                         << 
1065                                               << 
1066     for (G4int i = 1; i <= fNFiles; ++i) {  / << 
1067       rflag = std::fscanf(lecturePref, "%s",  << 
1068                                               << 
1069       G4String name = inputFile + G4String(". << 
1070       // Writes the results to a character st << 
1071                                               << 
1072       // G4cout << "check: " << name << G4end << 
1073       //  Open input file and give it to gest << 
1074       G4cout << "### Opening " << name << " a << 
1075       dicom = std::fopen(name.c_str(), "rb"); << 
1076       // Reading the .dcm in two steps:       << 
1077       //      1.  reading the header          << 
1078       //      2. reading the pixel data and s << 
1079       if (dicom != 0) {                       << 
1080         ReadFile(dicom, inputFile);           << 
1081       }                                       << 
1082       else {                                  << 
1083         G4cout << "\nError opening file : " <    944         G4cout << "\nError opening file : " << name << G4endl;
1084       }                                          945       }
1085       rflag = std::fclose(dicom);                946       rflag = std::fclose(dicom);
1086     }                                            947     }
1087 #endif                                        << 
1088                                               << 
1089     rflag = std::fclose(lecturePref);            948     rflag = std::fclose(lecturePref);
1090                                               << 949     
1091     // Checks the spacing is correct. Dumps t    950     // Checks the spacing is correct. Dumps to files
1092     fMergedSlices->CheckSlices();                951     fMergedSlices->CheckSlices();
1093                                               << 952     
1094     delete[] fCompressionc;                   << 953     delete [] fCompressionc;
1095     delete[] maxc;                            << 954     delete [] maxc;
1096     delete[] inputFile;                       << 955     delete [] name;
                                                   >> 956     delete [] inputFile;
1097     if (rflag) return;                           957     if (rflag) return;
                                                   >> 958     
1098   }                                              959   }
1099                                               << 960   
1100   if (fValueDensity) {                        << 961   if(fValueDensity) { delete [] fValueDensity; }
1101     delete[] fValueDensity;                   << 962   if(fValueCT) { delete [] fValueCT; }
1102   }                                           << 963   if(fMergedSlices) { delete fMergedSlices; }
1103   if (fValueCT) {                             << 964   
1104     delete[] fValueCT;                        << 
1105   }                                           << 
1106   if (fMergedSlices) {                        << 
1107     delete fMergedSlices;                     << 
1108   }                                           << 
1109 }                                                965 }
1110                                                  966 
1111 //....oooOO0OOooo........oooOO0OOooo........o << 967 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1112                                                  968 
1113 G4int DicomHandler::read_defined_nested(FILE* << 969 G4int DicomHandler::read_defined_nested(FILE * nested,G4int SQ_Length)
1114 {                                                970 {
1115   //      VARIABLES                              971   //      VARIABLES
1116   unsigned short item_GroupNumber;               972   unsigned short item_GroupNumber;
1117   unsigned short item_ElementNumber;             973   unsigned short item_ElementNumber;
1118   G4int item_Length;                             974   G4int item_Length;
1119   G4int items_array_length = 0;               << 975   G4int items_array_length=0;
1120   char* buffer = new char[LINEBUFFSIZE];      << 976   char * buffer= new char[LINEBUFFSIZE];
1121   size_t rflag = 0;                              977   size_t rflag = 0;
1122                                               << 978   
1123   while (items_array_length < SQ_Length) {    << 979   while(items_array_length < SQ_Length)
1124     rflag = std::fread(buffer, 2, 1, nested); << 980     {
1125     GetValue(buffer, item_GroupNumber);       << 981       rflag = std::fread(buffer, 2, 1, nested);
1126                                               << 982       GetValue(buffer, item_GroupNumber);
1127     rflag = std::fread(buffer, 2, 1, nested); << 983       
1128     GetValue(buffer, item_ElementNumber);     << 984       rflag = std::fread(buffer, 2, 1, nested);
1129                                               << 985       GetValue(buffer, item_ElementNumber);
1130     rflag = std::fread(buffer, 4, 1, nested); << 986       
1131     GetValue(buffer, item_Length);            << 987       rflag = std::fread(buffer, 4, 1, nested);
1132                                               << 988       GetValue(buffer, item_Length);
1133     rflag = std::fread(buffer, item_Length, 1 << 989       
1134                                               << 990       rflag = std::fread(buffer, item_Length, 1, nested);
1135     items_array_length = items_array_length + << 991       
1136   }                                           << 992       items_array_length= items_array_length+8+item_Length;
1137                                               << 993     }
1138   delete[] buffer;                            << 994   
1139                                               << 995   delete [] buffer;
1140   if (SQ_Length > items_array_length)         << 996   
                                                   >> 997   if( SQ_Length>items_array_length )
1141     return 0;                                    998     return 0;
1142   else                                           999   else
1143     return 1;                                    1000     return 1;
1144   if (rflag) return 1;                           1001   if (rflag) return 1;
1145 }                                                1002 }
1146                                                  1003 
1147 //....oooOO0OOooo........oooOO0OOooo........o << 1004 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1148                                                  1005 
1149 void DicomHandler::read_undefined_nested(FILE << 1006 void DicomHandler::read_undefined_nested(FILE * nested)
1150 {                                                1007 {
1151   //      VARIABLES                              1008   //      VARIABLES
1152   unsigned short item_GroupNumber;               1009   unsigned short item_GroupNumber;
1153   unsigned short item_ElementNumber;             1010   unsigned short item_ElementNumber;
1154   unsigned int item_Length;                      1011   unsigned int item_Length;
1155   char* buffer = new char[LINEBUFFSIZE];      << 1012   char * buffer= new char[LINEBUFFSIZE];
1156   size_t rflag = 0;                              1013   size_t rflag = 0;
1157                                               << 1014   
1158   do {                                        << 1015     do
1159     rflag = std::fread(buffer, 2, 1, nested); << 1016       {
1160     GetValue(buffer, item_GroupNumber);       << 1017         rflag = std::fread(buffer, 2, 1, nested);
1161                                               << 1018         GetValue(buffer, item_GroupNumber);
1162     rflag = std::fread(buffer, 2, 1, nested); << 1019         
1163     GetValue(buffer, item_ElementNumber);     << 1020         rflag = std::fread(buffer, 2, 1, nested);
1164                                               << 1021         GetValue(buffer, item_ElementNumber);
1165     rflag = std::fread(buffer, 4, 1, nested); << 1022         
1166     GetValue(buffer, item_Length);            << 1023         rflag = std::fread(buffer, 4, 1, nested);
1167                                               << 1024         GetValue(buffer, item_Length);
1168     if (item_Length != 0xffffffff)            << 1025         
1169       rflag = std::fread(buffer, item_Length, << 1026         if(item_Length!=0xffffffff)
1170     else                                      << 1027           rflag = std::fread(buffer, item_Length, 1, nested);
1171       read_undefined_item(nested);            << 1028         else
1172                                               << 1029           read_undefined_item(nested);
1173   } while (item_GroupNumber != 0xFFFE || item << 1030         
1174                                               << 1031         
1175   delete[] buffer;                            << 1032       } while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD 
1176   if (rflag) return;                          << 1033               || item_Length!=0);
                                                   >> 1034     
                                                   >> 1035     delete [] buffer;
                                                   >> 1036     if (rflag) return;
1177 }                                                1037 }
1178                                                  1038 
1179 //....oooOO0OOooo........oooOO0OOooo........o << 1039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1180                                                  1040 
1181 void DicomHandler::read_undefined_item(FILE*  << 1041 void DicomHandler::read_undefined_item(FILE * nested)
1182 {                                                1042 {
1183   //      VARIABLES                              1043   //      VARIABLES
1184   unsigned short item_GroupNumber;               1044   unsigned short item_GroupNumber;
1185   unsigned short item_ElementNumber;             1045   unsigned short item_ElementNumber;
1186   G4int item_Length;                          << 1046   G4int item_Length; size_t rflag = 0;
1187   size_t rflag = 0;                           << 1047   char *buffer= new char[LINEBUFFSIZE];
1188   char* buffer = new char[LINEBUFFSIZE];      << 1048   
1189                                               << 1049   do
1190   do {                                        << 1050     {
1191     rflag = std::fread(buffer, 2, 1, nested); << 1051       rflag = std::fread(buffer, 2, 1, nested);
1192     GetValue(buffer, item_GroupNumber);       << 1052       GetValue(buffer, item_GroupNumber);
1193                                               << 1053       
1194     rflag = std::fread(buffer, 2, 1, nested); << 1054       rflag = std::fread(buffer, 2, 1, nested);
1195     GetValue(buffer, item_ElementNumber);     << 1055       GetValue(buffer, item_ElementNumber);
1196                                               << 1056       
1197     rflag = std::fread(buffer, 4, 1, nested); << 1057       rflag = std::fread(buffer, 4, 1, nested);
1198     GetValue(buffer, item_Length);            << 1058       GetValue(buffer, item_Length);
1199                                               << 1059       
1200     if (item_Length != 0) rflag = std::fread( << 1060       
1201                                               << 1061       if(item_Length!=0)
1202   } while (item_GroupNumber != 0xFFFE || item << 1062         rflag = std::fread(buffer,item_Length,1,nested);
1203                                               << 1063       
1204   delete[] buffer;                            << 1064     }
                                                   >> 1065   while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D 
                                                   >> 1066         || item_Length!=0);
                                                   >> 1067   
                                                   >> 1068   delete [] buffer;
1205   if (rflag) return;                             1069   if (rflag) return;
1206 }                                                1070 }
1207                                                  1071 
1208 //....oooOO0OOooo........oooOO0OOooo........o << 1072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1209                                                  1073 
1210 template<class Type>                          << 1074 template <class Type>
1211 void DicomHandler::GetValue(char* _val, Type& << 1075 void DicomHandler::GetValue(char * _val, Type & _rval) {
1212 {                                             << 1076   
1213 #if BYTE_ORDER == BIG_ENDIAN                     1077 #if BYTE_ORDER == BIG_ENDIAN
1214   if (fLittleEndian) {  // little endian      << 1078   if(fLittleEndian) {      // little endian
1215 #else  // BYTE_ORDER == LITTLE_ENDIAN         << 1079 #else // BYTE_ORDER == LITTLE_ENDIAN
1216   if (!fLittleEndian) {  // big endian        << 1080     if(!fLittleEndian) {     // big endian
1217 #endif                                           1081 #endif
1218     const G4int SIZE = sizeof(_rval);         << 1082       const int SIZE = sizeof(_rval);
1219     char ctemp;                               << 1083       char ctemp;
1220     for (G4int i = 0; i < SIZE / 2; ++i) {    << 1084       for(int i = 0; i < SIZE/2; i++) {
1221       ctemp = _val[i];                        << 1085         ctemp = _val[i];
1222       _val[i] = _val[SIZE - 1 - i];           << 1086         _val[i] = _val[SIZE - 1 - i];
1223       _val[SIZE - 1 - i] = ctemp;             << 1087         _val[SIZE - 1 - i] = ctemp;
                                                   >> 1088       }
1224     }                                            1089     }
                                                   >> 1090     _rval = *(Type *)_val;
1225   }                                              1091   }
1226   _rval = *(Type*)_val;                       << 1092   
1227 }                                             << 1093   //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1228                                                  1094 
1229 //....oooOO0OOooo........oooOO0OOooo........o << 
1230                                                  1095