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 9.0)


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