| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file medical/DICOM/src/DicomHandler.cc
28 /// \brief Implementation of the DicomHandler class
29 //
30 // The code was written by :
31 // *Louis Archambault louis.archambault@phy.ulaval.ca,
32 // *Luc Beaulieu beaulieu@phy.ulaval.ca
33 // +Vincent Hubert-Tremblay at tigre.2@sympatico.ca
34 //
35 //
36 // *Centre Hospitalier Universitaire de Quebec (CHUQ),
37 // Hotel-Dieu de Quebec, departement de Radio-oncologie
38 // 11 cote du palais. Quebec, QC, Canada, G1R 2J6
39 // tel (418) 525-4444 #6720
40 // fax (418) 691 5268
41 //
42 // + University Laval, Quebec (QC) Canada
43 //*******************************************************
44 //
45 //*******************************************************
46 //
47 /// DicomHandler.cc :
48 /// - Handling of DICM images
49 /// - Reading headers and pixels
50 /// - Transforming pixel to density and creating *.g4dcm
51 /// files
52 //*******************************************************
53
54 #include "DicomHandler.hh"
55
56 #include "DicomPhantomZSliceHeader.hh"
57 #include "DicomPhantomZSliceMerged.hh"
58
59 #include "G4ios.hh"
60 #include "globals.hh"
61
62 #include <cctype>
63 #include <cstring>
64 #include <fstream>
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67
68 DicomHandler* DicomHandler::fInstance = 0;
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72 DicomHandler* DicomHandler::Instance()
73 {
74 if (fInstance == 0) {
75 static DicomHandler dicomhandler;
76 fInstance = &dicomhandler;
77 }
78 return fInstance;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........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 installation
96 #ifdef DICOM_USE_HEAD
97 G4cerr << "Warning! DICOM was compiled with DICOM_USE_HEAD option but "
98 << "the DICOM_PATH was not set!" << G4endl;
99 G4String datadir = G4GetEnv<G4String>("G4ENSDFSTATEDATA", "");
100 if (datadir.length() > 0) {
101 auto _last = datadir.rfind("/");
102 if (_last != std::string::npos) datadir.erase(_last);
103 driverPath = datadir + "/DICOM1.1/DICOM_HEAD_TEST";
104 G4int rc = setenv("DICOM_PATH", driverPath.c_str(), 0);
105 G4cerr << "\t --> Using '" << driverPath << "'..." << G4endl;
106 G4ConsumeParameters(rc);
107 }
108 #endif
109 }
110 return driverPath;
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114
115 G4String DicomHandler::GetDicomDataFile()
116 {
117 #if defined(DICOM_USE_HEAD) && defined(G4_DCMTK)
118 return GetDicomDataPath() + "/Data.dat.new";
119 #else
120 return GetDicomDataPath() + "/Data.dat";
121 #endif
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 " << fDriverFile << G4endl;
157 }
158 #else
159 DicomHandler::DicomHandler()
160 : DATABUFFSIZE(8192),
161 LINEBUFFSIZE(5020),
162 FILENAMESIZE(512),
163 fCompression(0),
164 fNFiles(0),
165 fRows(0),
166 fColumns(0),
167 fBitAllocated(0),
168 fMaxPixelValue(0),
169 fMinPixelValue(0),
170 fPixelSpacingX(0.),
171 fPixelSpacingY(0.),
172 fSliceThickness(0.),
173 fSliceLocation(0.),
174 fRescaleIntercept(0),
175 fRescaleSlope(0),
176 fLittleEndian(true),
177 fImplicitEndian(false),
178 fPixelRepresentation(0),
179 fNbrequali(0),
180 fValueDensity(NULL),
181 fValueCT(NULL),
182 fReadCalibration(false),
183 fMergedSlices(NULL),
184 fDriverFile("Data.dat"),
185 fCt2DensityFile("CT2Density.dat")
186 {
187 fMergedSlices = new DicomPhantomZSliceMerged;
188 }
189 #endif
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192 DicomHandler::~DicomHandler() {}
193
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
195
196 G4int DicomHandler::ReadFile(FILE* dicom, char* filename2)
197 {
198 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
207 rflag = std::fread(buffer, 1, 128, dicom); // The first 128 bytes
208 // are not important
209 // Reads the "DICOM" letters
210 rflag = std::fread(buffer, 1, 4, dicom);
211 // if there is no preamble, the FILE pointer is rewinded.
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 input data
218 short readElementId; // identify a particular type information
219 short elementLength2; // deal with element length in 2 bytes
220 // unsigned int elementLength4;
221 // deal with element length in 4 bytes
222 unsigned long elementLength4; // deal with element length in 4 bytes
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 afterward
239 G4int tagDictionary = readGroupId * 0x10000 + readElementId;
240
241 // beginning of the pixels
242 if (tagDictionary == 0x7FE00010) {
243 // Following 2 fread's are modifications to original DICOM example
244 // (Jonathan Madsen)
245 if (!fImplicitEndian) rflag = std::fread(buffer, 2, 1, dicom); // Reserved 2 bytes
246 // (not used for pixels)
247 rflag = std::fread(buffer, 4, 1, dicom); // Element Length
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, SQ, UN, added OF and UT
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); // Skip 2 reserved bytes
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, elementLength4) == 0) {
281 G4Exception("DicomHandler::ReadFile", "DICOM001", FatalException,
282 "Function read_defined_nested() failed!");
283 }
284 }
285 }
286 else {
287 // Reading the information with data
288 rflag = std::fread(data, elementLength4, 1, dicom);
289 }
290 }
291 else {
292 // explicit with VR different than previous ones
293 if (!fImplicitEndian || readGroupId == 2) {
294 // G4cout << "Reading DICOM files with Explicit VR"<< G4endl;
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, elementLength4, 1, dicom);
301 }
302 else { // Implicit VR
303
304 // G4cout << "Reading DICOM files with Implicit VR"<< G4endl;
305
306 // element length (4 bytes)
307 if (std::fseek(dicom, -2, SEEK_CUR) != 0) {
308 G4Exception("DicomHandler::ReadFile", "DICOM001", FatalException, "fseek failed");
309 }
310
311 rflag = std::fread(buffer, 4, 1, dicom);
312 GetValue(buffer, elementLength4);
313
314 // G4cout << std::hex<< elementLength4 << G4endl;
315
316 if (elementLength4 == 0xFFFFFFFF) {
317 read_undefined_nested(dicom);
318 elementLength4 = 0;
319 }
320 else {
321 rflag = std::fread(data, elementLength4, 1, dicom);
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) + ".g4dcm";
334
335 // Perform functions originally written straight to file
336 DicomPhantomZSliceHeader* zslice = new DicomPhantomZSliceHeader(fnameG4DCM);
337 for (auto ite = fMaterialIndices.cbegin(); ite != fMaterialIndices.cend(); ++ite) {
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 / 2.);
346 zslice->SetMaxX(fPixelSpacingX * fColumns / 2.);
347
348 zslice->SetMinY(-fPixelSpacingY * fRows / 2.);
349 zslice->SetMaxY(fPixelSpacingY * fRows / 2.);
350
351 zslice->SetMinZ(fSliceLocation - fSliceThickness / 2.);
352 zslice->SetMaxZ(fSliceLocation + fSliceThickness / 2.);
353
354 //===
355
356 ReadData(dicom, filename2);
357
358 // DEPRECIATED
359 // StoreData( foutG4DCM );
360 // foutG4DCM.close();
361
362 StoreData(zslice);
363
364 // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency
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 }
376
377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378
379 void DicomHandler::GetInformation(G4int& tagDictionary, char* data)
380 {
381 if (tagDictionary == 0x00280010) { // Number of Rows
382 GetValue(data, fRows);
383 std::printf("[0x00280010] Rows -> %i\n", fRows);
384 }
385 else if (tagDictionary == 0x00280011) { // Number of fColumns
386 GetValue(data, fColumns);
387 std::printf("[0x00280011] Columns -> %i\n", fColumns);
388 }
389 else if (tagDictionary == 0x00280102) { // High bits ( not used )
390 short highBits;
391 GetValue(data, highBits);
392 std::printf("[0x00280102] High bits -> %i\n", highBits);
393 }
394 else if (tagDictionary == 0x00280100) { // Bits allocated
395 GetValue(data, fBitAllocated);
396 std::printf("[0x00280100] Bits allocated -> %i\n", fBitAllocated);
397 }
398 else if (tagDictionary == 0x00280101) { // Bits stored ( not used )
399 short bitStored;
400 GetValue(data, bitStored);
401 std::printf("[0x00280101] Bits stored -> %i\n", bitStored);
402 }
403 else if (tagDictionary == 0x00280106) { // Min. pixel value
404 GetValue(data, fMinPixelValue);
405 std::printf("[0x00280106] Min. pixel value -> %i\n", fMinPixelValue);
406 }
407 else if (tagDictionary == 0x00280107) { // Max. pixel value
408 GetValue(data, fMaxPixelValue);
409 std::printf("[0x00280107] Max. pixel value -> %i\n", fMaxPixelValue);
410 }
411 else if (tagDictionary == 0x00281053) { // Rescale slope
412 fRescaleSlope = atoi(data);
413 std::printf("[0x00281053] Rescale Slope -> %d\n", fRescaleSlope);
414 }
415 else if (tagDictionary == 0x00281052) { // Rescalse intercept
416 fRescaleIntercept = atoi(data);
417 std::printf("[0x00281052] Rescale Intercept -> %d\n", fRescaleIntercept);
418 }
419 else if (tagDictionary == 0x00280103) {
420 // Pixel representation ( functions not design to read signed bits )
421 fPixelRepresentation = atoi(data); // 0: unsigned 1: signed
422 std::printf("[0x00280103] Pixel Representation -> %i\n", fPixelRepresentation);
423 if (fPixelRepresentation == 1) {
424 std::printf("### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
425 std::printf("DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
426 std::printf("ERROR !!!!!! -> \n");
427 }
428 }
429 else if (tagDictionary == 0x00080006) { // Modality
430 std::printf("[0x00080006] Modality -> %s\n", data);
431 }
432 else if (tagDictionary == 0x00080070) { // Manufacturer
433 std::printf("[0x00080070] Manufacturer -> %s\n", data);
434 }
435 else if (tagDictionary == 0x00080080) { // Institution Name
436 std::printf("[0x00080080] Institution Name -> %s\n", data);
437 }
438 else if (tagDictionary == 0x00080081) { // Institution Address
439 std::printf("[0x00080081] Institution Address -> %s\n", data);
440 }
441 else if (tagDictionary == 0x00081040) { // Institution Department Name
442 std::printf("[0x00081040] Institution Department Name -> %s\n", data);
443 }
444 else if (tagDictionary == 0x00081090) { // Manufacturer's Model Name
445 std::printf("[0x00081090] Manufacturer's Model Name -> %s\n", data);
446 }
447 else if (tagDictionary == 0x00181000) { // Device Serial Number
448 std::printf("[0x00181000] Device Serial Number -> %s\n", data);
449 }
450 else if (tagDictionary == 0x00080008) { // Image type ( not used )
451 std::printf("[0x00080008] Image Types -> %s\n", data);
452 }
453 else if (tagDictionary == 0x00283000) { // Modality LUT Sequence(not used)
454 std::printf("[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
455 }
456 else if (tagDictionary == 0x00283002) { // LUT Descriptor ( not used )
457 std::printf("[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
458 }
459 else if (tagDictionary == 0x00283003) { // LUT Explanation ( not used )
460 std::printf("[0x00283003] LUT Explanation LO 1 -> %s\n", data);
461 }
462 else if (tagDictionary == 0x00283004) { // Modality LUT ( not used )
463 std::printf("[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
464 }
465 else if (tagDictionary == 0x00283006) { // LUT Data ( not used )
466 std::printf("[0x00283006] LUT Data US or SS -> %s\n", data);
467 }
468 else if (tagDictionary == 0x00283010) { // VOI LUT ( not used )
469 std::printf("[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
470 }
471 else if (tagDictionary == 0x00280120) { // Pixel Padding Value (not used)
472 std::printf("[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
473 }
474 else if (tagDictionary == 0x00280030) { // Pixel Spacing
475 G4String datas(data);
476 G4int iss = G4int(datas.find('\\'));
477 fPixelSpacingX = atof(datas.substr(0, iss).c_str());
478 fPixelSpacingY = atof(datas.substr(iss + 1, datas.length()).c_str());
479 }
480 else if (tagDictionary == 0x00200037) { // Image Orientation ( not used )
481 std::printf("[0x00200037] Image Orientation (Phantom) -> %s\n", data);
482 }
483 else if (tagDictionary == 0x00200032) { // Image Position ( not used )
484 std::printf("[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
485 }
486 else if (tagDictionary == 0x00180050) { // Slice Thickness
487 fSliceThickness = atof(data);
488 std::printf("[0x00180050] Slice Thickness (mm) -> %f\n", fSliceThickness);
489 }
490 else if (tagDictionary == 0x00201041) { // Slice Location
491 fSliceLocation = atof(data);
492 std::printf("[0x00201041] Slice Location -> %f\n", fSliceLocation);
493 }
494 else if (tagDictionary == 0x00280004) { // Photometric Interpretation
495 // ( not used )
496 std::printf("[0x00280004] Photometric Interpretation -> %s\n", data);
497 }
498 else if (tagDictionary == 0x00020010) { // Endian
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.2", 19) == 0)
502 fLittleEndian = false;
503 // else 1.2.840..10008.1.2.1 (explicit little endian)
504
505 std::printf("[0x00020010] Endian -> %s\n", data);
506 }
507
508 // others
509 else {
510 // std::printf("[0x%x] -> %s\n", tagDictionary, data);
511 ;
512 }
513 }
514
515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
516
517 void DicomHandler::StoreData(DicomPhantomZSliceHeader* dcmPZSH)
518 {
519 G4int mean;
520 G4double density;
521 G4bool overflow = false;
522
523 if (!dcmPZSH) {
524 return;
525 }
526
527 dcmPZSH->SetSliceLocation(fSliceLocation);
528
529 //----- Print indices of material
530 if (fCompression == 1) { // no fCompression: each pixel has a density value)
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(G4float(density)));
538 }
539 }
540 }
541 else {
542 // density value is the average of a square region of
543 // fCompression*fCompression pixels
544 for (G4int ww = 0; ww < fRows; ww += fCompression) {
545 dcmPZSH->AddRow();
546 for (G4int xx = 0; xx < fColumns; xx += fCompression) {
547 overflow = false;
548 mean = 0;
549 for (G4int sumx = 0; sumx < fCompression; ++sumx) {
550 for (G4int sumy = 0; sumy < fCompression; ++sumy) {
551 if (ww + sumy >= fRows || xx + sumx >= fColumns) overflow = true;
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(G4float(density)));
562 }
563 }
564 }
565 }
566
567 dcmPZSH->FlipData();
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
571 // This function is depreciated as it is handled by
572 // DicomPhantomZSliceHeader::DumpToFile
573 void DicomHandler::StoreData(std::ofstream& foutG4DCM)
574 {
575 G4int mean;
576 G4double density;
577 G4bool overflow = false;
578
579 //----- Print indices of material
580 if (fCompression == 1) { // no fCompression: each pixel has a density value)
581 for (G4int ww = 0; ww < fRows; ++ww) {
582 for (G4int xx = 0; xx < fColumns; ++xx) {
583 mean = fTab[ww][xx];
584 density = Pixel2density(mean);
585 foutG4DCM << GetMaterialIndex(G4float(density)) << " ";
586 }
587 foutG4DCM << G4endl;
588 }
589 }
590 else {
591 // density value is the average of a square region of
592 // fCompression*fCompression pixels
593 for (G4int ww = 0; ww < fRows; ww += fCompression) {
594 for (G4int xx = 0; xx < fColumns; xx += fCompression) {
595 overflow = false;
596 mean = 0;
597 for (G4int sumx = 0; sumx < fCompression; ++sumx) {
598 for (G4int sumy = 0; sumy < fCompression; ++sumy) {
599 if (ww + sumy >= fRows || xx + sumx >= fColumns) overflow = true;
600 mean += fTab[ww + sumy][xx + sumx];
601 }
602 if (overflow) break;
603 }
604 mean /= fCompression * fCompression;
605
606 if (!overflow) {
607 density = Pixel2density(mean);
608 foutG4DCM << GetMaterialIndex(G4float(density)) << " ";
609 }
610 }
611 foutG4DCM << G4endl;
612 }
613 }
614
615 //----- Print densities
616 if (fCompression == 1) { // no fCompression: each pixel has a density value)
617 for (G4int ww = 0; ww < fRows; ww++) {
618 for (G4int xx = 0; xx < fColumns; xx++) {
619 mean = fTab[ww][xx];
620 density = Pixel2density(mean);
621 foutG4DCM << density << " ";
622 if (xx % 8 == 3) foutG4DCM << G4endl; // just for nicer reading
623 }
624 }
625 }
626 else {
627 // density value is the average of a square region of
628 // fCompression*fCompression pixels
629 for (G4int ww = 0; ww < fRows; ww += fCompression) {
630 for (G4int xx = 0; xx < fColumns; xx += fCompression) {
631 overflow = false;
632 mean = 0;
633 for (G4int sumx = 0; sumx < fCompression; ++sumx) {
634 for (G4int sumy = 0; sumy < fCompression; ++sumy) {
635 if (ww + sumy >= fRows || xx + sumx >= fColumns) overflow = true;
636 mean += fTab[ww + sumy][xx + sumx];
637 }
638 if (overflow) break;
639 }
640 mean /= fCompression * fCompression;
641
642 if (!overflow) {
643 density = Pixel2density(mean);
644 foutG4DCM << density << " ";
645 if (xx / fCompression % 8 == 3) foutG4DCM << G4endl; // just for nicer
646 // reading
647 }
648 }
649 }
650 }
651 }
652
653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
654 void DicomHandler::ReadMaterialIndices(std::ifstream& finData)
655 {
656 unsigned int nMate;
657 G4String mateName;
658 G4float densityMax;
659 finData >> nMate;
660 if (finData.eof()) return;
661
662 G4cout << " ReadMaterialIndices " << nMate << G4endl;
663 for (unsigned int ii = 0; ii < nMate; ++ii) {
664 finData >> mateName >> densityMax;
665 fMaterialIndices[densityMax] = mateName;
666 // G4cout << ii << " ReadMaterialIndices " << mateName << " "
667 // << densityMax << G4endl;
668 }
669 }
670
671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
672
673 unsigned int DicomHandler::GetMaterialIndex(G4float density)
674 {
675 std::map<G4float, G4String>::const_reverse_iterator ite;
676 std::size_t ii = fMaterialIndices.size();
677
678 for (ite = fMaterialIndices.crbegin(); ite != fMaterialIndices.crend(); ++ite, ii--) {
679 if (density >= (*ite).first) {
680 break;
681 }
682 }
683
684 if (ii == fMaterialIndices.size()) {
685 ii = fMaterialIndices.size() - 1;
686 }
687
688 return unsigned(ii);
689 }
690
691 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
692
693 G4int DicomHandler::ReadData(FILE* dicom, char* filename2)
694 {
695 G4int returnvalue = 0;
696 size_t rflag = 0;
697
698 // READING THE PIXELS
699
700 fTab = new G4int*[fRows];
701 for (G4int i = 0; i < fRows; ++i) {
702 fTab[i] = new G4int[fColumns];
703 }
704
705 if (fBitAllocated == 8) { // Case 8 bits :
706
707 std::printf("@@@ Error! Picture != 16 bits...\n");
708 std::printf("@@@ Error! Picture != 16 bits...\n");
709 std::printf("@@@ Error! Picture != 16 bits...\n");
710
711 unsigned char ch = 0;
712
713 for (G4int j = 0; j < fRows; ++j) {
714 for (G4int i = 0; i < fColumns; ++i) {
715 rflag = std::fread(&ch, 1, 1, dicom);
716 fTab[j][i] = ch * fRescaleSlope + fRescaleIntercept;
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 + fRescaleIntercept;
729 }
730 }
731 }
732
733 // Creation of .g4 files wich contains averaged density data
734 G4String nameProcessed = filename2 + G4String(".g4dcmb");
735 FILE* fileOut = std::fopen(nameProcessed.c_str(), "w+b");
736
737 G4cout << "### Writing of " << nameProcessed << " ###\n";
738
739 unsigned int nMate = fMaterialIndices.size();
740 rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
741 //--- Write materials
742 for (auto ite = fMaterialIndices.cbegin(); ite != fMaterialIndices.cend(); ++ite) {
743 G4String mateName = (*ite).second;
744 for (std::size_t ii = (*ite).second.length(); ii < 40; ++ii) {
745 mateName += " ";
746 } // mateName = const_cast<char*>(((*ite).second).c_str());
747
748 const char* mateNameC = mateName.c_str();
749 rflag = std::fwrite(mateNameC, sizeof(char), 40, fileOut);
750 }
751
752 unsigned int fRowsC = fRows / fCompression;
753 unsigned int fColumnsC = fColumns / fCompression;
754 unsigned int planesC = 1;
755 G4float pixelLocationXM = -G4float(fPixelSpacingX * fColumns / 2.);
756 G4float pixelLocationXP = G4float(fPixelSpacingX * fColumns / 2.);
757 G4float pixelLocationYM = -G4float(fPixelSpacingY * fRows / 2.);
758 G4float pixelLocationYP = G4float(fPixelSpacingY * fRows / 2.);
759 G4float fSliceLocationZM = G4float(fSliceLocation - fSliceThickness / 2.);
760 G4float fSliceLocationZP = G4float(fSliceLocation + fSliceThickness / 2.);
761 //--- Write number of voxels (assume only one voxel in Z)
762 rflag = std::fwrite(&fRowsC, sizeof(unsigned int), 1, fileOut);
763 rflag = std::fwrite(&fColumnsC, sizeof(unsigned int), 1, fileOut);
764 rflag = std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut);
765 //--- Write minimum and maximum extensions
766 rflag = std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut);
767 rflag = std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut);
768 rflag = std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut);
769 rflag = std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut);
770 rflag = std::fwrite(&fSliceLocationZM, sizeof(G4float), 1, fileOut);
771 rflag = std::fwrite(&fSliceLocationZP, sizeof(G4float), 1, fileOut);
772 // rflag = std::fwrite(&fCompression, sizeof(unsigned int), 1, fileOut);
773
774 std::printf("%8i %8i\n", fRows, fColumns);
775 std::printf("%8f %8f\n", fPixelSpacingX, fPixelSpacingY);
776 std::printf("%8f\n", fSliceThickness);
777 std::printf("%8f\n", fSliceLocation);
778 std::printf("%8i\n", fCompression);
779
780 G4int compSize = fCompression;
781 G4int mean;
782 G4float density;
783 G4bool overflow = false;
784
785 //----- Write index of material for each pixel
786 if (compSize == 1) { // no fCompression: each pixel has a density value)
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(density);
792 rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
793 }
794 }
795 }
796 else {
797 // density value is the average of a square region of
798 // fCompression*fCompression pixels
799 for (G4int ww = 0; ww < fRows; ww += compSize) {
800 for (G4int xx = 0; xx < fColumns; xx += compSize) {
801 overflow = false;
802 mean = 0;
803 for (G4int sumx = 0; sumx < compSize; ++sumx) {
804 for (G4int sumy = 0; sumy < compSize; ++sumy) {
805 if (ww + sumy >= fRows || xx + sumx >= fColumns) overflow = true;
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 = GetMaterialIndex(density);
815 rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
816 }
817 }
818 }
819 }
820
821 //----- Write density for each pixel
822 if (compSize == 1) { // no fCompression: each pixel has a density value)
823 for (G4int ww = 0; ww < fRows; ++ww) {
824 for (G4int xx = 0; xx < fColumns; ++xx) {
825 mean = fTab[ww][xx];
826 density = Pixel2density(mean);
827 rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
828 }
829 }
830 }
831 else {
832 // density value is the average of a square region of
833 // fCompression*fCompression pixels
834 for (G4int ww = 0; ww < fRows; ww += compSize) {
835 for (G4int xx = 0; xx < fColumns; xx += compSize) {
836 overflow = false;
837 mean = 0;
838 for (G4int sumx = 0; sumx < compSize; ++sumx) {
839 for (G4int sumy = 0; sumy < compSize; ++sumy) {
840 if (ww + sumy >= fRows || xx + sumx >= fColumns) overflow = true;
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(G4float), 1, fileOut);
850 }
851 }
852 }
853 }
854
855 rflag = std::fclose(fileOut);
856
857 delete[] nameProcessed;
858
859 /* for ( G4int i = 0; i < fRows; i ++ ) {
860 delete [] fTab[i];
861 }
862 delete [] fTab;
863 */
864
865 if (rflag) return returnvalue;
866 return returnvalue;
867 }
868
869 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......
870
871 // DICOM HEAD does not use the calibration curve
872
873 #ifdef DICOM_USE_HEAD
874 void DicomHandler::ReadCalibration()
875 {
876 fNbrequali = 0;
877 fReadCalibration = false;
878 G4cout << "No calibration curve for the DICOM HEAD needed!" << G4endl;
879 }
880 #else
881 // Separated out of Pixel2density
882 // No need to read in same calibration EVERY time
883 // Increases the speed of reading file by several orders of magnitude
884
885 void DicomHandler::ReadCalibration()
886 {
887 fNbrequali = 0;
888 // CT2Density.dat contains the calibration curve to convert CT (Hounsfield)
889 // number to physical density
890 std::ifstream calibration(fCt2DensityFile.c_str());
891 calibration >> fNbrequali;
892 fValueDensity = new G4double[fNbrequali];
893 fValueCT = new G4double[fNbrequali];
894
895 if (!calibration) {
896 G4Exception("DicomHandler::ReadFile", "DICOM001", FatalException,
897 "@@@ No value to transform pixels in density!");
898 }
899 else { // calibration was successfully opened
900 for (G4int i = 0; i < fNbrequali; ++i) { // Loop to store all the pts in CT2Density.dat
901 calibration >> fValueCT[i] >> fValueDensity[i];
902 }
903 }
904 calibration.close();
905 fReadCalibration = true;
906 }
907 #endif
908
909 #ifdef DICOM_USE_HEAD
910 G4float DicomHandler::Pixel2density(G4int pixel)
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........oooOO0OOooo.....
943
944 G4float DicomHandler::Pixel2density(G4int pixel)
945 {
946 if (!fReadCalibration) {
947 ReadCalibration();
948 }
949
950 G4float density = -1.;
951 G4double deltaCT = 0;
952 G4double deltaDensity = 0;
953
954 for (G4int j = 1; j < fNbrequali; ++j) {
955 if (pixel >= fValueCT[j - 1] && pixel < fValueCT[j]) {
956 deltaCT = fValueCT[j] - fValueCT[j - 1];
957 deltaDensity = fValueDensity[j] - fValueDensity[j - 1];
958
959 // interpolating linearly
960 density = G4float(fValueDensity[j] - ((fValueCT[j] - pixel) * deltaDensity / deltaCT));
961 break;
962 }
963 }
964
965 if (density < 0.) {
966 std::printf(
967 "@@@ Error density = %f && Pixel = %i \
968 (0x%x) && deltaDensity/deltaCT = %f\n",
969 density, pixel, pixel, deltaDensity / deltaCT);
970 }
971
972 return density;
973 }
974 #endif
975 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
976
977 void DicomHandler::CheckFileFormat()
978 {
979 std::ifstream checkData(fDriverFile.c_str());
980 char* oneLine = new char[128];
981
982 if (!(checkData.is_open())) { // Check existance of Data.dat
983
984 G4String message = "\nDicomG4 needs Data.dat (or another driver file specified";
985 message += " in command line):\n";
986 message += "\tFirst line: number of image pixel for a voxel (G4Box)\n";
987 message += "\tSecond line: number of images (CT slices) to read\n";
988 message += "\tEach following line contains the name of a Dicom image";
989 message += " except for the .dcm extension";
990 G4Exception("DicomHandler::ReadFile", "DICOM001", FatalException, message.c_str());
991 }
992
993 checkData >> fCompression;
994 checkData >> fNFiles;
995 G4String oneName;
996 checkData.getline(oneLine, 100);
997 std::ifstream testExistence;
998 G4bool existAlready = true;
999 for (G4int rep = 0; rep < fNFiles; ++rep) {
1000 checkData.getline(oneLine, 100);
1001 oneName = oneLine;
1002 oneName += ".g4dcm"; // create dicomFile.g4dcm
1003 G4cout << fNFiles << " test file " << oneName << G4endl;
1004 testExistence.open(oneName.data());
1005 if (!(testExistence.is_open())) {
1006 existAlready = false;
1007 testExistence.clear();
1008 testExistence.close();
1009 }
1010 testExistence.clear();
1011 testExistence.close();
1012 }
1013
1014 ReadMaterialIndices(checkData);
1015
1016 checkData.close();
1017 delete[] oneLine;
1018
1019 if (existAlready == false) { // The files *.g4dcm have to be created
1020
1021 G4cout << "\nAll the necessary images were not found in processed form "
1022 << ", starting with .dcm images\n";
1023
1024 FILE* dicom;
1025 FILE* lecturePref;
1026 char* fCompressionc = new char[LINEBUFFSIZE];
1027 char* maxc = new char[LINEBUFFSIZE];
1028 // char name[300], inputFile[300];
1029 char* inputFile = new char[FILENAMESIZE];
1030 G4int rflag;
1031 lecturePref = std::fopen(fDriverFile.c_str(), "r");
1032
1033 rflag = std::fscanf(lecturePref, "%s", fCompressionc);
1034 fCompression = atoi(fCompressionc);
1035 rflag = std::fscanf(lecturePref, "%s", maxc);
1036 fNFiles = atoi(maxc);
1037 G4cout << " fNFiles " << fNFiles << G4endl;
1038
1039 /////////////////////
1040
1041 #ifdef DICOM_USE_HEAD
1042 for (G4int i = 1; i <= fNFiles; ++i) { // Begin loop on filenames
1043 rflag = std::fscanf(lecturePref, "%s", inputFile);
1044 G4String path = GetDicomDataPath() + "/";
1045
1046 G4String name = inputFile + G4String(".dcm");
1047 // Writes the results to a character string buffer.
1048
1049 G4String name2 = path + name;
1050 // Open input file and give it to gestion_dicom :
1051 G4cout << "### Opening " << name2 << " and reading :\n";
1052 dicom = std::fopen(name2.c_str(), "rb");
1053 // Reading the .dcm in two steps:
1054 // 1. reading the header
1055 // 2. reading the pixel data and store the density in Moyenne.dat
1056 if (dicom != 0) {
1057 ReadFile(dicom, inputFile);
1058 }
1059 else {
1060 G4cout << "\nError opening file : " << name2 << G4endl;
1061 }
1062 rflag = std::fclose(dicom);
1063 }
1064 #else
1065
1066 for (G4int i = 1; i <= fNFiles; ++i) { // Begin loop on filenames
1067 rflag = std::fscanf(lecturePref, "%s", inputFile);
1068
1069 G4String name = inputFile + G4String(".dcm");
1070 // Writes the results to a character string buffer.
1071
1072 // G4cout << "check: " << name << G4endl;
1073 // Open input file and give it to gestion_dicom :
1074 G4cout << "### Opening " << name << " and reading :\n";
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 store the density in Moyenne.dat
1079 if (dicom != 0) {
1080 ReadFile(dicom, inputFile);
1081 }
1082 else {
1083 G4cout << "\nError opening file : " << name << G4endl;
1084 }
1085 rflag = std::fclose(dicom);
1086 }
1087 #endif
1088
1089 rflag = std::fclose(lecturePref);
1090
1091 // Checks the spacing is correct. Dumps to files
1092 fMergedSlices->CheckSlices();
1093
1094 delete[] fCompressionc;
1095 delete[] maxc;
1096 delete[] inputFile;
1097 if (rflag) return;
1098 }
1099
1100 if (fValueDensity) {
1101 delete[] fValueDensity;
1102 }
1103 if (fValueCT) {
1104 delete[] fValueCT;
1105 }
1106 if (fMergedSlices) {
1107 delete fMergedSlices;
1108 }
1109 }
1110
1111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1112
1113 G4int DicomHandler::read_defined_nested(FILE* nested, G4int SQ_Length)
1114 {
1115 // VARIABLES
1116 unsigned short item_GroupNumber;
1117 unsigned short item_ElementNumber;
1118 G4int item_Length;
1119 G4int items_array_length = 0;
1120 char* buffer = new char[LINEBUFFSIZE];
1121 size_t rflag = 0;
1122
1123 while (items_array_length < SQ_Length) {
1124 rflag = std::fread(buffer, 2, 1, nested);
1125 GetValue(buffer, item_GroupNumber);
1126
1127 rflag = std::fread(buffer, 2, 1, nested);
1128 GetValue(buffer, item_ElementNumber);
1129
1130 rflag = std::fread(buffer, 4, 1, nested);
1131 GetValue(buffer, item_Length);
1132
1133 rflag = std::fread(buffer, item_Length, 1, nested);
1134
1135 items_array_length = items_array_length + 8 + item_Length;
1136 }
1137
1138 delete[] buffer;
1139
1140 if (SQ_Length > items_array_length)
1141 return 0;
1142 else
1143 return 1;
1144 if (rflag) return 1;
1145 }
1146
1147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1148
1149 void DicomHandler::read_undefined_nested(FILE* nested)
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
1168 if (item_Length != 0xffffffff)
1169 rflag = std::fread(buffer, item_Length, 1, nested);
1170 else
1171 read_undefined_item(nested);
1172
1173 } while (item_GroupNumber != 0xFFFE || item_ElementNumber != 0xE0DD || item_Length != 0);
1174
1175 delete[] buffer;
1176 if (rflag) return;
1177 }
1178
1179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1180
1181 void DicomHandler::read_undefined_item(FILE* nested)
1182 {
1183 // VARIABLES
1184 unsigned short item_GroupNumber;
1185 unsigned short item_ElementNumber;
1186 G4int item_Length;
1187 size_t rflag = 0;
1188 char* buffer = new char[LINEBUFFSIZE];
1189
1190 do {
1191 rflag = std::fread(buffer, 2, 1, nested);
1192 GetValue(buffer, item_GroupNumber);
1193
1194 rflag = std::fread(buffer, 2, 1, nested);
1195 GetValue(buffer, item_ElementNumber);
1196
1197 rflag = std::fread(buffer, 4, 1, nested);
1198 GetValue(buffer, item_Length);
1199
1200 if (item_Length != 0) rflag = std::fread(buffer, item_Length, 1, nested);
1201
1202 } while (item_GroupNumber != 0xFFFE || item_ElementNumber != 0xE00D || item_Length != 0);
1203
1204 delete[] buffer;
1205 if (rflag) return;
1206 }
1207
1208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1209
1210 template<class Type>
1211 void DicomHandler::GetValue(char* _val, Type& _rval)
1212 {
1213 #if BYTE_ORDER == BIG_ENDIAN
1214 if (fLittleEndian) { // little endian
1215 #else // BYTE_ORDER == LITTLE_ENDIAN
1216 if (!fLittleEndian) { // big endian
1217 #endif
1218 const G4int SIZE = sizeof(_rval);
1219 char ctemp;
1220 for (G4int i = 0; i < SIZE / 2; ++i) {
1221 ctemp = _val[i];
1222 _val[i] = _val[SIZE - 1 - i];
1223 _val[SIZE - 1 - i] = ctemp;
1224 }
1225 }
1226 _rval = *(Type*)_val;
1227 }
1228
1229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1230