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