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