| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file medical/DICOM/src/DicomHandler.cc
28 /// \brief Implementation of the DicomHandler
29 //
30 // The code was written by :
31 // *Louis Archambault louis.archambault@p
32 // *Luc Beaulieu beaulieu@phy.ulaval.ca
33 // +Vincent Hubert-Tremblay at tigre.2@sy
34 //
35 //
36 // *Centre Hospitalier Universitaire de Quebec
37 // Hotel-Dieu de Quebec, departement de Radio-
38 // 11 cote du palais. Quebec, QC, Canada, G1R
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
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........oo
67
68 DicomHandler* DicomHandler::fInstance = 0;
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 DicomHandler* DicomHandler::Instance()
73 {
74 if (fInstance == 0) {
75 static DicomHandler dicomhandler;
76 fInstance = &dicomhandler;
77 }
78 return fInstance;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4String DicomHandler::GetDicomDataPath()
84 {
85 // default is current directory
86 G4String driverPath = ".";
87 // check environment
88 const char* path = std::getenv("DICOM_PATH")
89
90 if (path) {
91 // if is set in environment
92 return G4String(path);
93 }
94 else {
95 // if DICOM_USE_HEAD, look for data instal
96 #ifdef DICOM_USE_HEAD
97 G4cerr << "Warning! DICOM was compiled wit
98 << "the DICOM_PATH was not set!" <<
99 G4String datadir = G4GetEnv<G4String>("G4E
100 if (datadir.length() > 0) {
101 auto _last = datadir.rfind("/");
102 if (_last != std::string::npos) datadir.
103 driverPath = datadir + "/DICOM1.1/DICOM_
104 G4int rc = setenv("DICOM_PATH", driverPa
105 G4cerr << "\t --> Using '" << driverPath
106 G4ConsumeParameters(rc);
107 }
108 #endif
109 }
110 return driverPath;
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 G4String DicomHandler::GetDicomDataFile()
116 {
117 #if defined(DICOM_USE_HEAD) && defined(G4_DCMT
118 return GetDicomDataPath() + "/Data.dat.new";
119 #else
120 return GetDicomDataPath() + "/Data.dat";
121 #endif
122 }
123
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()
160 : DATABUFFSIZE(8192),
161 LINEBUFFSIZE(5020),
162 FILENAMESIZE(512),
163 fCompression(0),
164 fNFiles(0),
165 fRows(0),
166 fColumns(0),
167 fBitAllocated(0),
168 fMaxPixelValue(0),
169 fMinPixelValue(0),
170 fPixelSpacingX(0.),
171 fPixelSpacingY(0.),
172 fSliceThickness(0.),
173 fSliceLocation(0.),
174 fRescaleIntercept(0),
175 fRescaleSlope(0),
176 fLittleEndian(true),
177 fImplicitEndian(false),
178 fPixelRepresentation(0),
179 fNbrequali(0),
180 fValueDensity(NULL),
181 fValueCT(NULL),
182 fReadCalibration(false),
183 fMergedSlices(NULL),
184 fDriverFile("Data.dat"),
185 fCt2DensityFile("CT2Density.dat")
186 {
187 fMergedSlices = new DicomPhantomZSliceMerged
188 }
189 #endif
190 //....oooOO0OOooo........oooOO0OOooo........oo
191
192 DicomHandler::~DicomHandler() {}
193
194 //....oooOO0OOooo........oooOO0OOooo........oo
195
196 G4int DicomHandler::ReadFile(FILE* dicom, char
197 {
198 G4cout << " ReadFile " << filename2 << G4end
199
200 G4int returnvalue = 0;
201 size_t rflag = 0;
202 char* buffer = new char[LINEBUFFSIZE];
203
204 fImplicitEndian = false;
205 fLittleEndian = true;
206
207 rflag = std::fread(buffer, 1, 128, dicom);
208
209
210 rflag = std::fread(buffer, 1, 4, dicom);
211 // if there is no preamble, the FILE pointer
212 if (std::strncmp("DICM", buffer, 4) != 0) {
213 std::fseek(dicom, 0, SEEK_SET);
214 fImplicitEndian = true;
215 }
216
217 short readGroupId; // identify the kind of
218 short readElementId; // identify a particul
219 short elementLength2; // deal with element
220 // unsigned int eleme
221 // deal with element length in 4 bytes
222 unsigned long elementLength4; // deal with
223
224 char* data = new char[DATABUFFSIZE];
225
226 // Read information up to the pixel data
227 while (true) {
228 // Reading groups and elements :
229 readGroupId = 0;
230 readElementId = 0;
231 // group ID
232 rflag = std::fread(buffer, 2, 1, dicom);
233 GetValue(buffer, readGroupId);
234 // element ID
235 rflag = std::fread(buffer, 2, 1, dicom);
236 GetValue(buffer, readElementId);
237
238 // Creating a tag to be identified afterwa
239 G4int tagDictionary = readGroupId * 0x1000
240
241 // beginning of the pixels
242 if (tagDictionary == 0x7FE00010) {
243 // Following 2 fread's are modifications
244 // (Jonathan Madsen)
245 if (!fImplicitEndian) rflag = std::fread
246 // (not used for pixels)
247 rflag = std::fread(buffer, 4, 1, dicom);
248 // (not used for pixels)
249 break; // Exit to ReadImageData()
250 }
251
252 // VR or element length
253 rflag = std::fread(buffer, 2, 1, dicom);
254 GetValue(buffer, elementLength2);
255
256 // If value representation (VR) is OB, OW,
257 // the next length is 32 bits
258 if ((elementLength2 == 0x424f || // "OB"
259 elementLength2 == 0x574f || // "OW"
260 elementLength2 == 0x464f || // "OF"
261 elementLength2 == 0x5455 || // "UT"
262 elementLength2 == 0x5153 || // "SQ"
263 elementLength2 == 0x4e55)
264 && // "UN"
265 !fImplicitEndian)
266 { // explicit VR
267
268 rflag = std::fread(buffer, 2, 1, dicom);
269
270 // element length
271 rflag = std::fread(buffer, 4, 1, dicom);
272 GetValue(buffer, elementLength4);
273
274 if (elementLength2 == 0x5153) {
275 if (elementLength4 == 0xFFFFFFFF) {
276 read_undefined_nested(dicom);
277 elementLength4 = 0;
278 }
279 else {
280 if (read_defined_nested(dicom, eleme
281 G4Exception("DicomHandler::ReadFil
282 "Function read_defined
283 }
284 }
285 }
286 else {
287 // Reading the information with data
288 rflag = std::fread(data, elementLength
289 }
290 }
291 else {
292 // explicit with VR different than prev
293 if (!fImplicitEndian || readGroupId == 2
294 // G4cout << "Reading DICOM files wit
295 // element length (2 bytes)
296 rflag = std::fread(buffer, 2, 1, dicom
297 GetValue(buffer, elementLength2);
298 elementLength4 = elementLength2;
299
300 rflag = std::fread(data, elementLength
301 }
302 else { // Implicit VR
303
304 // G4cout << "Reading DICOM files wit
305
306 // element length (4 bytes)
307 if (std::fseek(dicom, -2, SEEK_CUR) !=
308 G4Exception("DicomHandler::ReadFile"
309 }
310
311 rflag = std::fread(buffer, 4, 1, dicom
312 GetValue(buffer, elementLength4);
313
314 // G4cout << std::hex<< elementLength
315
316 if (elementLength4 == 0xFFFFFFFF) {
317 read_undefined_nested(dicom);
318 elementLength4 = 0;
319 }
320 else {
321 rflag = std::fread(data, elementLeng
322 }
323 }
324 }
325
326 // NULL termination
327 data[elementLength4] = '\0';
328
329 // analyzing information
330 GetInformation(tagDictionary, data);
331 }
332
333 G4String fnameG4DCM = G4String(filename2) +
334
335 // Perform functions originally written stra
336 DicomPhantomZSliceHeader* zslice = new Dicom
337 for (auto ite = fMaterialIndices.cbegin(); i
338 zslice->AddMaterial(ite->second);
339 }
340
341 zslice->SetNoVoxelsX(fColumns / fCompression
342 zslice->SetNoVoxelsY(fRows / fCompression);
343 zslice->SetNoVoxelsZ(1);
344
345 zslice->SetMinX(-fPixelSpacingX * fColumns /
346 zslice->SetMaxX(fPixelSpacingX * fColumns /
347
348 zslice->SetMinY(-fPixelSpacingY * fRows / 2.
349 zslice->SetMaxY(fPixelSpacingY * fRows / 2.)
350
351 zslice->SetMinZ(fSliceLocation - fSliceThick
352 zslice->SetMaxZ(fSliceLocation + fSliceThick
353
354 //===
355
356 ReadData(dicom, filename2);
357
358 // DEPRECIATED
359 // StoreData( foutG4DCM );
360 // foutG4DCM.close();
361
362 StoreData(zslice);
363
364 // Dumped 2 file after DicomPhantomZSliceMer
365 // zslice->DumpToFile();
366
367 fMergedSlices->AddZSlice(zslice);
368
369 //
370 delete[] buffer;
371 delete[] data;
372
373 if (rflag) return returnvalue;
374 return returnvalue;
375 }
376
377 //....oooOO0OOooo........oooOO0OOooo........oo
378
379 void DicomHandler::GetInformation(G4int& tagDi
380 {
381 if (tagDictionary == 0x00280010) { // Numbe
382 GetValue(data, fRows);
383 std::printf("[0x00280010] Rows -> %i\n", f
384 }
385 else if (tagDictionary == 0x00280011) { //
386 GetValue(data, fColumns);
387 std::printf("[0x00280011] Columns -> %i\n"
388 }
389 else if (tagDictionary == 0x00280102) { //
390 short highBits;
391 GetValue(data, highBits);
392 std::printf("[0x00280102] High bits -> %i\
393 }
394 else if (tagDictionary == 0x00280100) { //
395 GetValue(data, fBitAllocated);
396 std::printf("[0x00280100] Bits allocated -
397 }
398 else if (tagDictionary == 0x00280101) { //
399 short bitStored;
400 GetValue(data, bitStored);
401 std::printf("[0x00280101] Bits stored -> %
402 }
403 else if (tagDictionary == 0x00280106) { //
404 GetValue(data, fMinPixelValue);
405 std::printf("[0x00280106] Min. pixel value
406 }
407 else if (tagDictionary == 0x00280107) { //
408 GetValue(data, fMaxPixelValue);
409 std::printf("[0x00280107] Max. pixel value
410 }
411 else if (tagDictionary == 0x00281053) { //
412 fRescaleSlope = atoi(data);
413 std::printf("[0x00281053] Rescale Slope ->
414 }
415 else if (tagDictionary == 0x00281052) { //
416 fRescaleIntercept = atoi(data);
417 std::printf("[0x00281052] Rescale Intercep
418 }
419 else if (tagDictionary == 0x00280103) {
420 // Pixel representation ( functions not d
421 fPixelRepresentation = atoi(data); // 0:
422 std::printf("[0x00280103] Pixel Representa
423 if (fPixelRepresentation == 1) {
424 std::printf("### PIXEL REPRESENTATION =
425 std::printf("DICOM READING SCAN FOR UNSI
426 std::printf("ERROR !!!!!! -> \n");
427 }
428 }
429 else if (tagDictionary == 0x00080006) { //
430 std::printf("[0x00080006] Modality -> %s\n
431 }
432 else if (tagDictionary == 0x00080070) { //
433 std::printf("[0x00080070] Manufacturer ->
434 }
435 else if (tagDictionary == 0x00080080) { //
436 std::printf("[0x00080080] Institution Name
437 }
438 else if (tagDictionary == 0x00080081) { //
439 std::printf("[0x00080081] Institution Addr
440 }
441 else if (tagDictionary == 0x00081040) { //
442 std::printf("[0x00081040] Institution Depa
443 }
444 else if (tagDictionary == 0x00081090) { //
445 std::printf("[0x00081090] Manufacturer's M
446 }
447 else if (tagDictionary == 0x00181000) { //
448 std::printf("[0x00181000] Device Serial Nu
449 }
450 else if (tagDictionary == 0x00080008) { //
451 std::printf("[0x00080008] Image Types -> %
452 }
453 else if (tagDictionary == 0x00283000) { //
454 std::printf("[0x00283000] Modality LUT Seq
455 }
456 else if (tagDictionary == 0x00283002) { //
457 std::printf("[0x00283002] LUT Descriptor U
458 }
459 else if (tagDictionary == 0x00283003) { //
460 std::printf("[0x00283003] LUT Explanation
461 }
462 else if (tagDictionary == 0x00283004) { //
463 std::printf("[0x00283004] Modality LUT Typ
464 }
465 else if (tagDictionary == 0x00283006) { //
466 std::printf("[0x00283006] LUT Data US or S
467 }
468 else if (tagDictionary == 0x00283010) { //
469 std::printf("[0x00283010] VOI LUT Sequence
470 }
471 else if (tagDictionary == 0x00280120) { //
472 std::printf("[0x00280120] Pixel Padding Va
473 }
474 else if (tagDictionary == 0x00280030) { //
475 G4String datas(data);
476 G4int iss = G4int(datas.find('\\'));
477 fPixelSpacingX = atof(datas.substr(0, iss)
478 fPixelSpacingY = atof(datas.substr(iss + 1
479 }
480 else if (tagDictionary == 0x00200037) { //
481 std::printf("[0x00200037] Image Orientatio
482 }
483 else if (tagDictionary == 0x00200032) { //
484 std::printf("[0x00200032] Image Position (
485 }
486 else if (tagDictionary == 0x00180050) { //
487 fSliceThickness = atof(data);
488 std::printf("[0x00180050] Slice Thickness
489 }
490 else if (tagDictionary == 0x00201041) { //
491 fSliceLocation = atof(data);
492 std::printf("[0x00201041] Slice Location -
493 }
494 else if (tagDictionary == 0x00280004) { //
495 // ( not used )
496 std::printf("[0x00280004] Photometric Inte
497 }
498 else if (tagDictionary == 0x00020010) { //
499 if (strcmp(data, "1.2.840.10008.1.2") == 0
500 fImplicitEndian = true;
501 else if (strncmp(data, "1.2.840.10008.1.2.
502 fLittleEndian = false;
503 // else 1.2.840..10008.1.2.1 (explicit lit
504
505 std::printf("[0x00020010] Endian -> %s\n",
506 }
507
508 // others
509 else {
510 // std::printf("[0x%x] -> %s\n", tagDictio
511 ;
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
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 }
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
567 dcmPZSH->FlipData();
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oo
571 // This function is depreciated as it is handl
572 // DicomPhantomZSliceHeader::DumpToFile
573 void DicomHandler::StoreData(std::ofstream& fo
574 {
575 G4int mean;
576 G4double density;
577 G4bool overflow = false;
578
579 //----- Print indices of material
580 if (fCompression == 1) { // no fCompression
581 for (G4int ww = 0; ww < fRows; ++ww) {
582 for (G4int xx = 0; xx < fColumns; ++xx)
583 mean = fTab[ww][xx];
584 density = Pixel2density(mean);
585 foutG4DCM << GetMaterialIndex(G4float(
586 }
587 foutG4DCM << G4endl;
588 }
589 }
590 else {
591 // density value is the average of a squar
592 // fCompression*fCompression pixels
593 for (G4int ww = 0; ww < fRows; ww += fComp
594 for (G4int xx = 0; xx < fColumns; xx +=
595 overflow = false;
596 mean = 0;
597 for (G4int sumx = 0; sumx < fCompressi
598 for (G4int sumy = 0; sumy < fCompres
599 if (ww + sumy >= fRows || xx + sum
600 mean += fTab[ww + sumy][xx + sumx]
601 }
602 if (overflow) break;
603 }
604 mean /= fCompression * fCompression;
605
606 if (!overflow) {
607 density = Pixel2density(mean);
608 foutG4DCM << GetMaterialIndex(G4floa
609 }
610 }
611 foutG4DCM << G4endl;
612 }
613 }
614
615 //----- Print densities
616 if (fCompression == 1) { // no fCompression
617 for (G4int ww = 0; ww < fRows; ww++) {
618 for (G4int xx = 0; xx < fColumns; xx++)
619 mean = fTab[ww][xx];
620 density = Pixel2density(mean);
621 foutG4DCM << density << " ";
622 if (xx % 8 == 3) foutG4DCM << G4endl;
623 }
624 }
625 }
626 else {
627 // density value is the average of a squar
628 // fCompression*fCompression pixels
629 for (G4int ww = 0; ww < fRows; ww += fComp
630 for (G4int xx = 0; xx < fColumns; xx +=
631 overflow = false;
632 mean = 0;
633 for (G4int sumx = 0; sumx < fCompressi
634 for (G4int sumy = 0; sumy < fCompres
635 if (ww + sumy >= fRows || xx + sum
636 mean += fTab[ww + sumy][xx + sumx]
637 }
638 if (overflow) break;
639 }
640 mean /= fCompression * fCompression;
641
642 if (!overflow) {
643 density = Pixel2density(mean);
644 foutG4DCM << density << " ";
645 if (xx / fCompression % 8 == 3) fout
646 // reading
647 }
648 }
649 }
650 }
651 }
652
653 //....oooOO0OOooo........oooOO0OOooo........oo
654 void DicomHandler::ReadMaterialIndices(std::if
655 {
656 unsigned int nMate;
657 G4String mateName;
658 G4float densityMax;
659 finData >> nMate;
660 if (finData.eof()) return;
661
662 G4cout << " ReadMaterialIndices " << nMate <
663 for (unsigned int ii = 0; ii < nMate; ++ii)
664 finData >> mateName >> densityMax;
665 fMaterialIndices[densityMax] = mateName;
666 // G4cout << ii << " ReadMaterialIndice
667 // << densityMax << G4endl;
668 }
669 }
670
671 //....oooOO0OOooo........oooOO0OOooo........oo
672
673 unsigned int DicomHandler::GetMaterialIndex(G4
674 {
675 std::map<G4float, G4String>::const_reverse_i
676 std::size_t ii = fMaterialIndices.size();
677
678 for (ite = fMaterialIndices.crbegin(); ite !
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........oo
692
693 G4int DicomHandler::ReadData(FILE* dicom, char
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
708 std::printf("@@@ Error! Picture != 16 bits
709 std::printf("@@@ Error! Picture != 16 bits
710
711 unsigned char ch = 0;
712
713 for (G4int j = 0; j < fRows; ++j) {
714 for (G4int i = 0; i < fColumns; ++i) {
715 rflag = std::fread(&ch, 1, 1, dicom);
716 fTab[j][i] = ch * fRescaleSlope + fRes
717 }
718 }
719 returnvalue = 1;
720 }
721 else { // from 12 to 16 bits :
722 char sbuff[2];
723 short pixel;
724 for (G4int j = 0; j < fRows; ++j) {
725 for (G4int i = 0; i < fColumns; ++i) {
726 rflag = std::fread(sbuff, 2, 1, dicom)
727 GetValue(sbuff, pixel);
728 fTab[j][i] = pixel * fRescaleSlope + f
729 }
730 }
731 }
732
733 // Creation of .g4 files wich contains avera
734 G4String nameProcessed = filename2 + G4Strin
735 FILE* fileOut = std::fopen(nameProcessed.c_s
736
737 G4cout << "### Writing of " << nameProcessed
738
739 unsigned int nMate = fMaterialIndices.size()
740 rflag = std::fwrite(&nMate, sizeof(unsigned
741 //--- Write materials
742 for (auto ite = fMaterialIndices.cbegin(); i
743 G4String mateName = (*ite).second;
744 for (std::size_t ii = (*ite).second.length
745 mateName += " ";
746 } // mateName = const_cast<char*>(((*ite)
747
748 const char* mateNameC = mateName.c_str();
749 rflag = std::fwrite(mateNameC, sizeof(char
750 }
751
752 unsigned int fRowsC = fRows / fCompression;
753 unsigned int fColumnsC = fColumns / fCompres
754 unsigned int planesC = 1;
755 G4float pixelLocationXM = -G4float(fPixelSpa
756 G4float pixelLocationXP = G4float(fPixelSpac
757 G4float pixelLocationYM = -G4float(fPixelSpa
758 G4float pixelLocationYP = G4float(fPixelSpac
759 G4float fSliceLocationZM = G4float(fSliceLoc
760 G4float fSliceLocationZP = G4float(fSliceLoc
761 //--- Write number of voxels (assume only on
762 rflag = std::fwrite(&fRowsC, sizeof(unsigned
763 rflag = std::fwrite(&fColumnsC, sizeof(unsig
764 rflag = std::fwrite(&planesC, sizeof(unsigne
765 //--- Write minimum and maximum extensions
766 rflag = std::fwrite(&pixelLocationXM, sizeof
767 rflag = std::fwrite(&pixelLocationXP, sizeof
768 rflag = std::fwrite(&pixelLocationYM, sizeof
769 rflag = std::fwrite(&pixelLocationYP, sizeof
770 rflag = std::fwrite(&fSliceLocationZM, sizeo
771 rflag = std::fwrite(&fSliceLocationZP, sizeo
772 // rflag = std::fwrite(&fCompression, sizeof
773
774 std::printf("%8i %8i\n", fRows, fColumns);
775 std::printf("%8f %8f\n", fPixelSpacingX, f
776 std::printf("%8f\n", fSliceThickness);
777 std::printf("%8f\n", fSliceLocation);
778 std::printf("%8i\n", fCompression);
779
780 G4int compSize = fCompression;
781 G4int mean;
782 G4float density;
783 G4bool overflow = false;
784
785 //----- Write index of material for each pix
786 if (compSize == 1) { // no fCompression: ea
787 for (G4int ww = 0; ww < fRows; ++ww) {
788 for (G4int xx = 0; xx < fColumns; ++xx)
789 mean = fTab[ww][xx];
790 density = Pixel2density(mean);
791 unsigned int mateID = GetMaterialIndex
792 rflag = std::fwrite(&mateID, sizeof(un
793 }
794 }
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
821 //----- Write density for each pixel
822 if (compSize == 1) { // no fCompression: ea
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(G
828 }
829 }
830 }
831 else {
832 // density value is the average of a squar
833 // fCompression*fCompression pixels
834 for (G4int ww = 0; ww < fRows; ww += compS
835 for (G4int xx = 0; xx < fColumns; xx +=
836 overflow = false;
837 mean = 0;
838 for (G4int sumx = 0; sumx < compSize;
839 for (G4int sumy = 0; sumy < compSize
840 if (ww + sumy >= fRows || xx + sum
841 mean += fTab[ww + sumy][xx + sumx]
842 }
843 if (overflow) break;
844 }
845 mean /= compSize * compSize;
846
847 if (!overflow) {
848 density = Pixel2density(mean);
849 rflag = std::fwrite(&density, sizeof
850 }
851 }
852 }
853 }
854
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........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 }
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
950 G4float density = -1.;
951 G4double deltaCT = 0;
952 G4double deltaDensity = 0;
953
954 for (G4int j = 1; j < fNbrequali; ++j) {
955 if (pixel >= fValueCT[j - 1] && pixel < fV
956 deltaCT = fValueCT[j] - fValueCT[j - 1];
957 deltaDensity = fValueDensity[j] - fValue
958
959 // interpolating linearly
960 density = G4float(fValueDensity[j] - ((f
961 break;
962 }
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 / de
970 }
971
972 return density;
973 }
974 #endif
975 //....oooOO0OOooo........oooOO0OOooo........oo
976
977 void DicomHandler::CheckFileFormat()
978 {
979 std::ifstream checkData(fDriverFile.c_str())
980 char* oneLine = new char[128];
981
982 if (!(checkData.is_open())) { // Check exis
983
984 G4String message = "\nDicomG4 needs Data.d
985 message += " in command line):\n";
986 message += "\tFirst line: number of image
987 message += "\tSecond line: number of image
988 message += "\tEach following line contains
989 message += " except for the .dcm extension
990 G4Exception("DicomHandler::ReadFile", "DIC
991 }
992
993 checkData >> fCompression;
994 checkData >> fNFiles;
995 G4String oneName;
996 checkData.getline(oneLine, 100);
997 std::ifstream testExistence;
998 G4bool existAlready = true;
999 for (G4int rep = 0; rep < fNFiles; ++rep) {
1000 checkData.getline(oneLine, 100);
1001 oneName = oneLine;
1002 oneName += ".g4dcm"; // create dicomFile
1003 G4cout << fNFiles << " test file " << one
1004 testExistence.open(oneName.data());
1005 if (!(testExistence.is_open())) {
1006 existAlready = false;
1007 testExistence.clear();
1008 testExistence.close();
1009 }
1010 testExistence.clear();
1011 testExistence.close();
1012 }
1013
1014 ReadMaterialIndices(checkData);
1015
1016 checkData.close();
1017 delete[] oneLine;
1018
1019 if (existAlready == false) { // The files
1020
1021 G4cout << "\nAll the necessary images wer
1022 << ", starting with .dcm images\n"
1023
1024 FILE* dicom;
1025 FILE* lecturePref;
1026 char* fCompressionc = new char[LINEBUFFSI
1027 char* maxc = new char[LINEBUFFSIZE];
1028 // char name[300], inputFile[300];
1029 char* inputFile = new char[FILENAMESIZE];
1030 G4int rflag;
1031 lecturePref = std::fopen(fDriverFile.c_st
1032
1033 rflag = std::fscanf(lecturePref, "%s", fC
1034 fCompression = atoi(fCompressionc);
1035 rflag = std::fscanf(lecturePref, "%s", ma
1036 fNFiles = atoi(maxc);
1037 G4cout << " fNFiles " << fNFiles << G4end
1038
1039 /////////////////////
1040
1041 #ifdef DICOM_USE_HEAD
1042 for (G4int i = 1; i <= fNFiles; ++i) { /
1043 rflag = std::fscanf(lecturePref, "%s",
1044 G4String path = GetDicomDataPath() + "/
1045
1046 G4String name = inputFile + G4String(".
1047 // Writes the results to a character st
1048
1049 G4String name2 = path + name;
1050 // Open input file and give it to gest
1051 G4cout << "### Opening " << name2 << "
1052 dicom = std::fopen(name2.c_str(), "rb")
1053 // Reading the .dcm in two steps:
1054 // 1. reading the header
1055 // 2. reading the pixel data and s
1056 if (dicom != 0) {
1057 ReadFile(dicom, inputFile);
1058 }
1059 else {
1060 G4cout << "\nError opening file : " <
1061 }
1062 rflag = std::fclose(dicom);
1063 }
1064 #else
1065
1066 for (G4int i = 1; i <= fNFiles; ++i) { /
1067 rflag = std::fscanf(lecturePref, "%s",
1068
1069 G4String name = inputFile + G4String(".
1070 // Writes the results to a character st
1071
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
1091 // Checks the spacing is correct. Dumps t
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........o
1112
1113 G4int DicomHandler::read_defined_nested(FILE*
1114 {
1115 // VARIABLES
1116 unsigned short item_GroupNumber;
1117 unsigned short item_ElementNumber;
1118 G4int item_Length;
1119 G4int items_array_length = 0;
1120 char* buffer = new char[LINEBUFFSIZE];
1121 size_t rflag = 0;
1122
1123 while (items_array_length < SQ_Length) {
1124 rflag = std::fread(buffer, 2, 1, nested);
1125 GetValue(buffer, item_GroupNumber);
1126
1127 rflag = std::fread(buffer, 2, 1, nested);
1128 GetValue(buffer, item_ElementNumber);
1129
1130 rflag = std::fread(buffer, 4, 1, nested);
1131 GetValue(buffer, item_Length);
1132
1133 rflag = std::fread(buffer, item_Length, 1
1134
1135 items_array_length = items_array_length +
1136 }
1137
1138 delete[] buffer;
1139
1140 if (SQ_Length > items_array_length)
1141 return 0;
1142 else
1143 return 1;
1144 if (rflag) return 1;
1145 }
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
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 }
1178
1179 //....oooOO0OOooo........oooOO0OOooo........o
1180
1181 void DicomHandler::read_undefined_item(FILE*
1182 {
1183 // VARIABLES
1184 unsigned short item_GroupNumber;
1185 unsigned short item_ElementNumber;
1186 G4int item_Length;
1187 size_t rflag = 0;
1188 char* buffer = new char[LINEBUFFSIZE];
1189
1190 do {
1191 rflag = std::fread(buffer, 2, 1, nested);
1192 GetValue(buffer, item_GroupNumber);
1193
1194 rflag = std::fread(buffer, 2, 1, nested);
1195 GetValue(buffer, item_ElementNumber);
1196
1197 rflag = std::fread(buffer, 4, 1, nested);
1198 GetValue(buffer, item_Length);
1199
1200 if (item_Length != 0) rflag = std::fread(
1201
1202 } while (item_GroupNumber != 0xFFFE || item
1203
1204 delete[] buffer;
1205 if (rflag) return;
1206 }
1207
1208 //....oooOO0OOooo........oooOO0OOooo........o
1209
1210 template<class Type>
1211 void DicomHandler::GetValue(char* _val, Type&
1212 {
1213 #if BYTE_ORDER == BIG_ENDIAN
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........o
1230