Geant4  10.01.p03
DicomHandler.cc
Go to the documentation of this file.
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 // $Id: DicomHandler.cc 84839 2014-10-21 13:44:55Z gcosmo $
27 //
30 //
31 // The code was written by :
32 // *Louis Archambault louis.archambault@phy.ulaval.ca,
33 // *Luc Beaulieu beaulieu@phy.ulaval.ca
34 // +Vincent Hubert-Tremblay at tigre.2@sympatico.ca
35 //
36 //
37 // *Centre Hospitalier Universitaire de Quebec (CHUQ),
38 // Hotel-Dieu de Quebec, departement de Radio-oncologie
39 // 11 cote du palais. Quebec, QC, Canada, G1R 2J6
40 // tel (418) 525-4444 #6720
41 // fax (418) 691 5268
42 //
43 // + University Laval, Quebec (QC) Canada
44 //*******************************************************
45 //
46 //*******************************************************
47 //
53 //*******************************************************
54 
55 #include "DicomHandler.hh"
56 #include "globals.hh"
57 #include "G4ios.hh"
58 #include <fstream>
59 
60 #include <cctype>
61 #include <cstring>
62 
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
73 {
74  return fgInstance;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
80 : DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
81  fCompression(0), fNFiles(0), fRows(0), fColumns(0),
82  fBitAllocated(0), fMaxPixelValue(0), fMinPixelValue(0),
83  fPixelSpacingX(0.), fPixelSpacingY(0.),
84  fSliceThickness(0.), fSliceLocation(0.),
85  fRescaleIntercept(0), fRescaleSlope(0),
86  fLittleEndian(true), fImplicitEndian(false),
87  fPixelRepresentation(0), nbrequali(0),
88  valuedensity(NULL),valueCT(NULL),readCalibration(false),
89  mergedSlices(NULL),driverFile("Data.dat"),ct2densityFile("CT2Density.dat")
90 {
92 }
93 
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
98 {
99 
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 G4int DicomHandler::ReadFile(FILE* dicom, char* filename2)
105 {
106  G4cout << " ReadFile " << filename2 << G4endl;
107  G4int returnvalue = 0; size_t rflag = 0;
108  char * buffer = new char[LINEBUFFSIZE];
109 
110  fImplicitEndian = false;
111  fLittleEndian = true;
112 
113  rflag = std::fread( buffer, 1, 128, dicom ); // The first 128 bytes
114  //are not important
115  // Reads the "DICOM" letters
116  rflag = std::fread( buffer, 1, 4, dicom );
117  // if there is no preamble, the FILE pointer is rewinded.
118  if(std::strncmp("DICM", buffer, 4) != 0) {
119  std::fseek(dicom, 0, SEEK_SET);
120  fImplicitEndian = true;
121  }
122 
123  short readGroupId; // identify the kind of input data
124  short readElementId; // identify a particular type information
125  short elementLength2; // deal with element length in 2 bytes
126  //unsigned int elementLength4;
127  // deal with element length in 4 bytes
128  unsigned long elementLength4; // deal with element length in 4 bytes
129 
130  char * data = new char[DATABUFFSIZE];
131 
132 
133  // Read information up to the pixel data
134  while(true) {
135 
136  //Reading groups and elements :
137  readGroupId = 0;
138  readElementId = 0;
139  // group ID
140  rflag = std::fread(buffer, 2, 1, dicom);
141  GetValue(buffer, readGroupId);
142  // element ID
143  rflag = std::fread(buffer, 2, 1, dicom);
144  GetValue(buffer, readElementId);
145 
146  // Creating a tag to be identified afterward
147  G4int tagDictionary = readGroupId*0x10000 + readElementId;
148 
149  // beginning of the pixels
150  if(tagDictionary == 0x7FE00010) {
151  // Folling 2 fread's are modifications to original DICOM example
152  //(Jonathan Madsen)
153  rflag = std::fread(buffer,2,1,dicom); // Reserved 2 bytes
154  // (not used for pixels)
155  rflag = std::fread(buffer,4,1,dicom); // Element Length
156  // (not used for pixels)
157  break; // Exit to ReadImageData()
158  }
159 
160  // VR or element length
161  rflag = std::fread(buffer,2,1,dicom);
162  GetValue(buffer, elementLength2);
163 
164  // If value representation (VR) is OB, OW, SQ, UN, added OF and UT
165  //the next length is 32 bits
166  if((elementLength2 == 0x424f || // "OB"
167  elementLength2 == 0x574f || // "OW"
168  elementLength2 == 0x464f || // "OF"
169  elementLength2 == 0x5455 || // "UT"
170  elementLength2 == 0x5153 || // "SQ"
171  elementLength2 == 0x4e55) && // "UN"
172  !fImplicitEndian ) { // explicit VR
173 
174  rflag = std::fread(buffer, 2, 1, dicom); // Skip 2 reserved bytes
175 
176  // element length
177  rflag = std::fread(buffer, 4, 1, dicom);
178  GetValue(buffer, elementLength4);
179 
180  if(elementLength2 == 0x5153)
181  {
182  if(elementLength4 == 0xFFFFFFFF)
183  {
184  read_undefined_nested( dicom );
185  elementLength4=0;
186  } else{
187  if(read_defined_nested( dicom, elementLength4 )==0){
188  G4Exception("DicomHandler::ReadFile",
189  "DICOM001",
191  "Function read_defined_nested() failed!");
192  }
193  }
194  } else {
195  // Reading the information with data
196  rflag = std::fread(data, elementLength4,1,dicom);
197  }
198 
199 
200  } else {
201 
202  // explicit with VR different than previous ones
203  if(!fImplicitEndian || readGroupId == 2) {
204 
205  //G4cout << "Reading DICOM files with Explicit VR"<< G4endl;
206  // element length (2 bytes)
207  rflag = std::fread(buffer, 2, 1, dicom);
208  GetValue(buffer, elementLength2);
209  elementLength4 = elementLength2;
210 
211  rflag = std::fread(data, elementLength4, 1, dicom);
212 
213  } else { // Implicit VR
214 
215  //G4cout << "Reading DICOM files with Implicit VR"<< G4endl;
216 
217  // element length (4 bytes)
218  if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
219  G4Exception("DicomHandler::ReadFile",
220  "DICOM001",
222  "fseek failed");
223  }
224 
225  rflag = std::fread(buffer, 4, 1, dicom);
226  GetValue(buffer, elementLength4);
227 
228  //G4cout << std::hex<< elementLength4 << G4endl;
229 
230  if(elementLength4 == 0xFFFFFFFF)
231  {
232  read_undefined_nested(dicom);
233  elementLength4=0;
234  } else{
235  rflag = std::fread(data, elementLength4, 1, dicom);
236  }
237 
238  }
239  }
240 
241  // NULL termination
242  data[elementLength4] = '\0';
243 
244  // analyzing information
245  GetInformation(tagDictionary, data);
246  }
247 
248  G4String fnameG4DCM = G4String(filename2) + ".g4dcm";
249 
250  // Perform functions originally written straight to file
251  DicomPhantomZSliceHeader* zslice = new DicomPhantomZSliceHeader(fnameG4DCM);
252 
253  std::map<G4float,G4String>::const_iterator ite;
254  for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ++ite){
255  zslice->AddMaterial(ite->second);
256  }
257 
259  zslice->SetNoVoxelY(fRows/fCompression);
260  zslice->SetNoVoxelZ(1);
261 
262  zslice->SetMinX(-fPixelSpacingX*fColumns/2.);
263  zslice->SetMaxX(fPixelSpacingX*fColumns/2.);
264 
265  zslice->SetMinY(-fPixelSpacingY*fRows/2.);
266  zslice->SetMaxY(fPixelSpacingY*fRows/2.);
267 
270 
271  //===
272 
273  ReadData( dicom, filename2 );
274 
275  // DEPRECIATED
276  //StoreData( foutG4DCM );
277  //foutG4DCM.close();
278 
279  StoreData( zslice );
280 
281  // Dumped 2 file after DicomPhantomZSliceMerged has checked for consistency
282  //zslice->DumpToFile();
283 
284  mergedSlices->AddZSlice(zslice);
285 
286  //
287  delete [] buffer;
288  delete [] data;
289 
290  if (rflag) return returnvalue;
291  return returnvalue;
292 }
293 
294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
295 
296 void DicomHandler::GetInformation(G4int & tagDictionary, char * data)
297 {
298  if(tagDictionary == 0x00280010 ) { // Number of Rows
299  GetValue(data, fRows);
300  std::printf("[0x00280010] Rows -> %i\n",fRows);
301 
302  } else if(tagDictionary == 0x00280011 ) { // Number of fColumns
303  GetValue(data, fColumns);
304  std::printf("[0x00280011] Columns -> %i\n",fColumns);
305 
306  } else if(tagDictionary == 0x00280102 ) { // High bits ( not used )
307  short highBits;
308  GetValue(data, highBits);
309  std::printf("[0x00280102] High bits -> %i\n",highBits);
310 
311  } else if(tagDictionary == 0x00280100 ) { // Bits allocated
312  GetValue(data, fBitAllocated);
313  std::printf("[0x00280100] Bits allocated -> %i\n", fBitAllocated);
314 
315  } else if(tagDictionary == 0x00280101 ) { // Bits stored ( not used )
316  short bitStored;
317  GetValue(data, bitStored);
318  std::printf("[0x00280101] Bits stored -> %i\n",bitStored);
319 
320  } else if(tagDictionary == 0x00280106 ) { // Min. pixel value
321  GetValue(data, fMinPixelValue);
322  std::printf("[0x00280106] Min. pixel value -> %i\n", fMinPixelValue);
323 
324  } else if(tagDictionary == 0x00280107 ) { // Max. pixel value
325  GetValue(data, fMaxPixelValue);
326  std::printf("[0x00280107] Max. pixel value -> %i\n", fMaxPixelValue);
327 
328  } else if(tagDictionary == 0x00281053) { // Rescale slope
329  fRescaleSlope = atoi(data);
330  std::printf("[0x00281053] Rescale Slope -> %d\n", fRescaleSlope);
331 
332  } else if(tagDictionary == 0x00281052 ) { // Rescalse intercept
333  fRescaleIntercept = atoi(data);
334  std::printf("[0x00281052] Rescale Intercept -> %d\n",
336 
337  } else if(tagDictionary == 0x00280103 ) {
338  // Pixel representation ( functions not design to read signed bits )
339  fPixelRepresentation = atoi(data); // 0: unsigned 1: signed
340  std::printf("[0x00280103] Pixel Representation -> %i\n",
342  if(fPixelRepresentation == 1 ) {
343  std::printf("### PIXEL REPRESENTATION = 1, BITS ARE SIGNED, ");
344  std::printf("DICOM READING SCAN FOR UNSIGNED VALUE, POSSIBLE ");
345  std::printf("ERROR !!!!!! -> \n");
346  }
347 
348  } else if(tagDictionary == 0x00080006 ) { // Modality
349  std::printf("[0x00080006] Modality -> %s\n", data);
350 
351  } else if(tagDictionary == 0x00080070 ) { // Manufacturer
352  std::printf("[0x00080070] Manufacturer -> %s\n", data);
353 
354  } else if(tagDictionary == 0x00080080 ) { // Institution Name
355  std::printf("[0x00080080] Institution Name -> %s\n", data);
356 
357  } else if(tagDictionary == 0x00080081 ) { // Institution Address
358  std::printf("[0x00080081] Institution Address -> %s\n", data);
359 
360  } else if(tagDictionary == 0x00081040 ) { // Institution Department Name
361  std::printf("[0x00081040] Institution Department Name -> %s\n", data);
362 
363  } else if(tagDictionary == 0x00081090 ) { // Manufacturer's Model Name
364  std::printf("[0x00081090] Manufacturer's Model Name -> %s\n", data);
365 
366  } else if(tagDictionary == 0x00181000 ) { // Device Serial Number
367  std::printf("[0x00181000] Device Serial Number -> %s\n", data);
368 
369  } else if(tagDictionary == 0x00080008 ) { // Image type ( not used )
370  std::printf("[0x00080008] Image Types -> %s\n", data);
371 
372  } else if(tagDictionary == 0x00283000 ) { //Modality LUT Sequence(not used)
373  std::printf("[0x00283000] Modality LUT Sequence SQ 1 -> %s\n", data);
374 
375  } else if(tagDictionary == 0x00283002 ) { // LUT Descriptor ( not used )
376  std::printf("[0x00283002] LUT Descriptor US or SS 3 -> %s\n", data);
377 
378  } else if(tagDictionary == 0x00283003 ) { // LUT Explanation ( not used )
379  std::printf("[0x00283003] LUT Explanation LO 1 -> %s\n", data);
380 
381  } else if(tagDictionary == 0x00283004 ) { // Modality LUT ( not used )
382  std::printf("[0x00283004] Modality LUT Type LO 1 -> %s\n", data);
383 
384  } else if(tagDictionary == 0x00283006 ) { // LUT Data ( not used )
385  std::printf("[0x00283006] LUT Data US or SS -> %s\n", data);
386 
387  } else if(tagDictionary == 0x00283010 ) { // VOI LUT ( not used )
388  std::printf("[0x00283010] VOI LUT Sequence SQ 1 -> %s\n", data);
389 
390  } else if(tagDictionary == 0x00280120 ) { // Pixel Padding Value (not used)
391  std::printf("[0x00280120] Pixel Padding Value US or SS 1 -> %s\n", data);
392 
393  } else if(tagDictionary == 0x00280030 ) { // Pixel Spacing
394  G4String datas(data);
395  int iss = datas.find('\\');
396  fPixelSpacingX = atof( datas.substr(0,iss).c_str() );
397  fPixelSpacingY = atof( datas.substr(iss+2,datas.length()).c_str() );
398 
399  } else if(tagDictionary == 0x00200037 ) { // Image Orientation ( not used )
400  std::printf("[0x00200037] Image Orientation (Phantom) -> %s\n", data);
401 
402  } else if(tagDictionary == 0x00200032 ) { // Image Position ( not used )
403  std::printf("[0x00200032] Image Position (Phantom,mm) -> %s\n", data);
404 
405  } else if(tagDictionary == 0x00180050 ) { // Slice Thickness
406  fSliceThickness = atof(data);
407  std::printf("[0x00180050] Slice Thickness (mm) -> %f\n", fSliceThickness);
408 
409  } else if(tagDictionary == 0x00201041 ) { // Slice Location
410  fSliceLocation = atof(data);
411  std::printf("[0x00201041] Slice Location -> %f\n", fSliceLocation);
412 
413  } else if(tagDictionary == 0x00280004 ) { // Photometric Interpretation
414  // ( not used )
415  std::printf("[0x00280004] Photometric Interpretation -> %s\n", data);
416 
417  } else if(tagDictionary == 0x00020010) { // Endian
418  if(strcmp(data, "1.2.840.10008.1.2") == 0)
419  fImplicitEndian = true;
420  else if(strncmp(data, "1.2.840.10008.1.2.2", 19) == 0)
421  fLittleEndian = false;
422  //else 1.2.840..10008.1.2.1 (explicit little endian)
423 
424  std::printf("[0x00020010] Endian -> %s\n", data);
425  }
426 
427  // others
428  else {
429  //std::printf("[0x%x] -> %s\n", tagDictionary, data);
430  ;
431  }
432 
433 }
434 
435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
436 
438 {
439  G4int mean;
441  G4bool overflow = false;
442 
443  if(!dcmPZSH) { return; }
444 
446 
447  //----- Print indices of material
448  if(fCompression == 1) { // no fCompression: each pixel has a density value)
449  for( G4int ww = 0; ww < fRows; ww++) {
450  dcmPZSH->AddRow();
451  for( G4int xx = 0; xx < fColumns; xx++) {
452  mean = fTab[ww][xx];
453  density = Pixel2density(mean);
454  dcmPZSH->AddValue(density);
455  dcmPZSH->AddMateID(GetMaterialIndex(density));
456  }
457  }
458 
459  } else {
460  // density value is the average of a square region of
461  // fCompression*fCompression pixels
462  for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
463  dcmPZSH->AddRow();
464  for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
465  overflow = false;
466  mean = 0;
467  for(int sumx = 0; sumx < fCompression; sumx++) {
468  for(int sumy = 0; sumy < fCompression; sumy++) {
469  if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
470  mean += fTab[ww+sumy][xx+sumx];
471  }
472  if(overflow) break;
473  }
474  mean /= fCompression*fCompression;
475 
476  if(!overflow) {
477  density = Pixel2density(mean);
478  dcmPZSH->AddValue(density);
479  dcmPZSH->AddMateID(GetMaterialIndex(density));
480  }
481  }
482  }
483  }
484 
485  dcmPZSH->FlipData();
486 }
487 
488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489 // This function is depreciated as it is handled by
490 // DicomPhantomZSliceHeader::DumpToFile
491 void DicomHandler::StoreData(std::ofstream& foutG4DCM)
492 {
493  G4int mean;
495  G4bool overflow = false;
496 
497  //----- Print indices of material
498  if(fCompression == 1) { // no fCompression: each pixel has a density value)
499  for( G4int ww = 0; ww < fRows; ww++) {
500  for( G4int xx = 0; xx < fColumns; xx++) {
501  mean = fTab[ww][xx];
502  density = Pixel2density(mean);
503  foutG4DCM << GetMaterialIndex( density ) << " ";
504  }
505  foutG4DCM << G4endl;
506  }
507 
508  } else {
509  // density value is the average of a square region of
510  // fCompression*fCompression pixels
511  for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
512  for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
513  overflow = false;
514  mean = 0;
515  for(int sumx = 0; sumx < fCompression; sumx++) {
516  for(int sumy = 0; sumy < fCompression; sumy++) {
517  if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
518  mean += fTab[ww+sumy][xx+sumx];
519  }
520  if(overflow) break;
521  }
522  mean /= fCompression*fCompression;
523 
524  if(!overflow) {
525  density = Pixel2density(mean);
526  foutG4DCM << GetMaterialIndex( density ) << " ";
527  }
528  }
529  foutG4DCM << G4endl;
530  }
531 
532  }
533 
534  //----- Print densities
535  if(fCompression == 1) { // no fCompression: each pixel has a density value)
536  for( G4int ww = 0; ww < fRows; ww++) {
537  for( G4int xx = 0; xx < fColumns; xx++) {
538  mean = fTab[ww][xx];
539  density = Pixel2density(mean);
540  foutG4DCM << density << " ";
541  if( xx%8 == 3 ) foutG4DCM << G4endl; // just for nicer reading
542  }
543  }
544 
545  } else {
546  // density value is the average of a square region of
547  // fCompression*fCompression pixels
548  for(G4int ww = 0; ww < fRows ;ww += fCompression ) {
549  for(G4int xx = 0; xx < fColumns ;xx +=fCompression ) {
550  overflow = false;
551  mean = 0;
552  for(int sumx = 0; sumx < fCompression; sumx++) {
553  for(int sumy = 0; sumy < fCompression; sumy++) {
554  if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
555  mean += fTab[ww+sumy][xx+sumx];
556  }
557  if(overflow) break;
558  }
559  mean /= fCompression*fCompression;
560 
561  if(!overflow) {
562  density = Pixel2density(mean);
563  foutG4DCM << density << " ";
564  if( xx/fCompression%8 == 3 ) foutG4DCM << G4endl; // just for nicer
565  // reading
566  }
567  }
568  }
569 
570  }
571 
572 }
573 
574 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
575 void DicomHandler::ReadMaterialIndices( std::ifstream& finData)
576 {
577  unsigned int nMate;
578  G4String mateName;
579  G4float densityMax;
580  finData >> nMate;
581  if( finData.eof() ) return;
582 
583  G4cout << " ReadMaterialIndices " << nMate << G4endl;
584  for( unsigned int ii = 0; ii < nMate; ii++ ){
585  finData >> mateName >> densityMax;
586  fMaterialIndices[densityMax] = mateName;
587  G4cout << ii << " ReadMaterialIndices " << mateName << " "
588  << densityMax << G4endl;
589  }
590 
591 }
592 
593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
594 
596 {
597  std::map<G4float,G4String>::reverse_iterator ite;
598  G4int ii = fMaterialIndices.size();
599  for( ite = fMaterialIndices.rbegin(); ite != fMaterialIndices.rend();
600  ite++, ii-- ) {
601  if( density >= (*ite).first ) {
602  break;
603  }
604  }
605  //- G4cout << " GetMaterialIndex " << density << " = " << ii << G4endl;
606 
607  if(static_cast<unsigned int>(ii) == fMaterialIndices.size())
608  { ii = fMaterialIndices.size()-1; }
609 
610  return ii;
611 
612 }
613 
614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
615 
616 G4int DicomHandler::ReadData(FILE *dicom,char * filename2)
617 {
618  G4int returnvalue = 0; size_t rflag = 0;
619 
620  // READING THE PIXELS :
621  G4int w = 0;
622 
623  fTab = new G4int*[fRows];
624  for ( G4int i = 0; i < fRows; i ++ ) {
625  fTab[i] = new G4int[fColumns];
626  }
627 
628  if(fBitAllocated == 8) { // Case 8 bits :
629 
630  std::printf("@@@ Error! Picture != 16 bits...\n");
631  std::printf("@@@ Error! Picture != 16 bits...\n");
632  std::printf("@@@ Error! Picture != 16 bits...\n");
633 
634  unsigned char ch = 0;
635 
636  for(G4int j = 0; j < fRows; j++) {
637  for(G4int i = 0; i < fColumns; i++) {
638  w++;
639  rflag = std::fread( &ch, 1, 1, dicom);
640  fTab[j][i] = ch*fRescaleSlope + fRescaleIntercept;
641  }
642  }
643  returnvalue = 1;
644 
645  } else { // from 12 to 16 bits :
646  char sbuff[2];
647  short pixel;
648  for( G4int j = 0; j < fRows; j++) {
649  for( G4int i = 0; i < fColumns; i++) {
650  w++;
651  rflag = std::fread(sbuff, 2, 1, dicom);
652  GetValue(sbuff, pixel);
653  fTab[j][i] = pixel*fRescaleSlope + fRescaleIntercept;
654  }
655  }
656  }
657 
658  // Creation of .g4 files wich contains averaged density data
659  char * nameProcessed = new char[FILENAMESIZE];
660  FILE* fileOut;
661 
662  std::sprintf(nameProcessed,"%s.g4dcmb",filename2);
663  fileOut = std::fopen(nameProcessed,"w+b");
664  std::printf("### Writing of %s ###\n",nameProcessed);
665 
666  unsigned int nMate = fMaterialIndices.size();
667  rflag = std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
668  //--- Write materials
669  std::map<G4float,G4String>::const_iterator ite;
670  for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++ ){
671  G4String mateName = (*ite).second;
672  for( G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
673  mateName += " ";
674  } //mateName = const_cast<char*>(((*ite).second).c_str());
675 
676  const char* mateNameC = mateName.c_str();
677  rflag = std::fwrite(mateNameC, sizeof(char),40, fileOut);
678  }
679 
680  unsigned int fRowsC = fRows/fCompression;
681  unsigned int fColumnsC = fColumns/fCompression;
682  unsigned int planesC = 1;
683  G4float pixelLocationXM = -fPixelSpacingX*fColumns/2.;
684  G4float pixelLocationXP = fPixelSpacingX*fColumns/2.;
685  G4float pixelLocationYM = -fPixelSpacingY*fRows/2.;
686  G4float pixelLocationYP = fPixelSpacingY*fRows/2.;
687  G4float fSliceLocationZM = fSliceLocation-fSliceThickness/2.;
688  G4float fSliceLocationZP = fSliceLocation+fSliceThickness/2.;
689  //--- Write number of voxels (assume only one voxel in Z)
690  rflag = std::fwrite(&fRowsC, sizeof(unsigned int), 1, fileOut);
691  rflag = std::fwrite(&fColumnsC, sizeof(unsigned int), 1, fileOut);
692  rflag = std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut);
693  //--- Write minimum and maximum extensions
694  rflag = std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut);
695  rflag = std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut);
696  rflag = std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut);
697  rflag = std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut);
698  rflag = std::fwrite(&fSliceLocationZM, sizeof(G4float), 1, fileOut);
699  rflag = std::fwrite(&fSliceLocationZP, sizeof(G4float), 1, fileOut);
700  // rflag = std::fwrite(&fCompression, sizeof(unsigned int), 1, fileOut);
701 
702  std::printf("%8i %8i\n",fRows,fColumns);
703  std::printf("%8f %8f\n",fPixelSpacingX,fPixelSpacingY);
704  std::printf("%8f\n", fSliceThickness);
705  std::printf("%8f\n", fSliceLocation);
706  std::printf("%8i\n", fCompression);
707 
708  G4int compSize = fCompression;
709  G4int mean;
711  G4bool overflow = false;
712 
713  //----- Write index of material for each pixel
714  if(compSize == 1) { // no fCompression: each pixel has a density value)
715  for( G4int ww = 0; ww < fRows; ww++) {
716  for( G4int xx = 0; xx < fColumns; xx++) {
717  mean = fTab[ww][xx];
718  density = Pixel2density(mean);
719  unsigned int mateID = GetMaterialIndex( density );
720  rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
721  }
722  }
723 
724  } else {
725  // density value is the average of a square region of
726  // fCompression*fCompression pixels
727  for(G4int ww = 0; ww < fRows ;ww += compSize ) {
728  for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
729  overflow = false;
730  mean = 0;
731  for(int sumx = 0; sumx < compSize; sumx++) {
732  for(int sumy = 0; sumy < compSize; sumy++) {
733  if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
734  mean += fTab[ww+sumy][xx+sumx];
735  }
736  if(overflow) break;
737  }
738  mean /= compSize*compSize;
739 
740  if(!overflow) {
741  density = Pixel2density(mean);
742  unsigned int mateID = GetMaterialIndex( density );
743  rflag = std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
744  }
745  }
746 
747  }
748  }
749 
750  //----- Write density for each pixel
751  if(compSize == 1) { // no fCompression: each pixel has a density value)
752  for( G4int ww = 0; ww < fRows; ww++) {
753  for( G4int xx = 0; xx < fColumns; xx++) {
754  mean = fTab[ww][xx];
755  density = Pixel2density(mean);
756  rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
757  }
758  }
759 
760  } else {
761  // density value is the average of a square region of
762  // fCompression*fCompression pixels
763  for(G4int ww = 0; ww < fRows ;ww += compSize ) {
764  for(G4int xx = 0; xx < fColumns ;xx +=compSize ) {
765  overflow = false;
766  mean = 0;
767  for(int sumx = 0; sumx < compSize; sumx++) {
768  for(int sumy = 0; sumy < compSize; sumy++) {
769  if(ww+sumy >= fRows || xx+sumx >= fColumns) overflow = true;
770  mean += fTab[ww+sumy][xx+sumx];
771  }
772  if(overflow) break;
773  }
774  mean /= compSize*compSize;
775 
776  if(!overflow) {
777  density = Pixel2density(mean);
778  rflag = std::fwrite(&density, sizeof(G4float), 1, fileOut);
779  }
780  }
781 
782  }
783  }
784 
785  rflag = std::fclose(fileOut);
786 
787  delete [] nameProcessed;
788 
789  /* for ( G4int i = 0; i < fRows; i ++ ) {
790  delete [] fTab[i];
791  }
792  delete [] fTab;
793  */
794 
795  if (rflag) return returnvalue;
796  return returnvalue;
797 }
798 
799 
800 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
801 // Separated out of Pixel2density
802 // No need to read in same calibration EVERY time
803 // Increases the speed of reading file by several orders of magnitude
805 {
806  nbrequali = 0;
807 
808  // CT2Density.dat contains the calibration curve to convert CT (Hounsfield)
809  // number to physical density
810  std::ifstream calibration(ct2densityFile.c_str());
811  calibration >> nbrequali;
812 
814  valueCT = new G4double[nbrequali];
815 
816  if(!calibration) {
817  G4Exception("DicomHandler::ReadFile",
818  "DICOM001",
820  "@@@ No value to transform pixels in density!");
821 
822  } else { // calibration was successfully opened
823  for(G4int i = 0; i < nbrequali; i++) { // Loop to store all the
824  //pts in CT2Density.dat
825  calibration >> valueCT[i] >> valuedensity[i];
826  }
827  }
828  calibration.close();
829 
830  readCalibration = true;
831 }
832 
833 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
834 
836 {
837  if(!readCalibration) { ReadCalibration(); }
838 
839  G4float density = -1.;
840  G4double deltaCT = 0;
841  G4double deltaDensity = 0;
842 
843 
844  for(G4int j = 1; j < nbrequali; j++) {
845  if( pixel >= valueCT[j-1] && pixel < valueCT[j]) {
846 
847  deltaCT = valueCT[j] - valueCT[j-1];
848  deltaDensity = valuedensity[j] - valuedensity[j-1];
849 
850  // interpolating linearly
851  density = valuedensity[j] - ((valueCT[j] - pixel)*deltaDensity/deltaCT );
852  break;
853  }
854  }
855 
856  if(density < 0.) {
857  std::printf("@@@ Error density = %f && Pixel = %i \
858  (0x%x) && deltaDensity/deltaCT = %f\n",density,pixel,pixel,
859  deltaDensity/deltaCT);
860  }
861 
862  return density;
863 }
864 
865 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
866 
868 {
869  std::ifstream checkData(driverFile.c_str());
870  char * oneLine = new char[128];
871 
872  if(!(checkData.is_open())) { //Check existance of Data.dat
873 
874  G4String message =
875  "\nDicomG4 needs Data.dat (or another driver file specified";
876  message += " in command line):\n";
877  message += "\tFirst line: number of image pixel for a voxel (G4Box)\n";
878  message += "\tSecond line: number of images (CT slices) to read\n";
879  message += "\tEach following line contains the name of a Dicom image";
880  message += " except for the .dcm extension";
881  G4Exception("DicomHandler::ReadFile",
882  "DICOM001",
884  message.c_str());
885  }
886 
887  checkData >> fCompression;
888  checkData >> fNFiles;
889  G4String oneName;
890  checkData.getline(oneLine,100);
891  std::ifstream testExistence;
892  G4bool existAlready = true;
893  for(G4int rep = 0; rep < fNFiles; rep++) {
894  checkData.getline(oneLine,100);
895  oneName = oneLine;
896  oneName += ".g4dcm"; // create dicomFile.g4dcm
897  G4cout << fNFiles << " test file " << oneName << G4endl;
898  testExistence.open(oneName.data());
899  if(!(testExistence.is_open())) {
900  existAlready = false;
901  testExistence.clear();
902  testExistence.close();
903  }
904  testExistence.clear();
905  testExistence.close();
906  }
907 
908  ReadMaterialIndices( checkData );
909 
910  checkData.close();
911  delete [] oneLine;
912 
913  if( existAlready == false ) { // The files *.g4dcm have to be created
914 
915  G4cout << "\nAll the necessary images were not found in processed form "
916  << ", starting with .dcm images\n";
917 
918  FILE * dicom;
919  FILE * lecturePref;
920  char * fCompressionc = new char[LINEBUFFSIZE];
921  char * maxc = new char[LINEBUFFSIZE];
922  //char name[300], inputFile[300];
923  char * name = new char[FILENAMESIZE];
924  char * inputFile = new char[FILENAMESIZE];
925  G4int rflag;
926 
927  lecturePref = std::fopen(driverFile.c_str(),"r");
928  rflag = std::fscanf(lecturePref,"%s",fCompressionc);
929  fCompression = atoi(fCompressionc);
930  rflag = std::fscanf(lecturePref,"%s",maxc);
931  fNFiles = atoi(maxc);
932  G4cout << " fNFiles " << fNFiles << G4endl;
933 
934  for( G4int i = 1; i <= fNFiles; i++ ) { // Begin loop on filenames
935 
936  rflag = std::fscanf(lecturePref,"%s",inputFile);
937  std::sprintf(name,"%s.dcm",inputFile);
938  G4cout << "check 1: " << name << G4endl;
939  // Open input file and give it to gestion_dicom :
940  std::printf("### Opening %s and reading :\n",name);
941  dicom = std::fopen(name,"rb");
942  // Reading the .dcm in two steps:
943  // 1. reading the header
944  // 2. reading the pixel data and store the density in Moyenne.dat
945  if( dicom != 0 ) {
946  ReadFile(dicom,inputFile);
947  } else {
948  G4cout << "\nError opening file : " << name << G4endl;
949  }
950  rflag = std::fclose(dicom);
951  }
952  rflag = std::fclose(lecturePref);
953 
954  // Checks the spacing is correct. Dumps to files
956 
957  delete [] fCompressionc;
958  delete [] maxc;
959  delete [] name;
960  delete [] inputFile;
961  if (rflag) return;
962 
963  }
964 
965  if(valuedensity) { delete [] valuedensity; }
966  if(valueCT) { delete [] valueCT; }
967  if(mergedSlices) { delete mergedSlices; }
968 
969 }
970 
971 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
972 
974 {
975  // VARIABLES
976  unsigned short item_GroupNumber;
977  unsigned short item_ElementNumber;
978  G4int item_Length;
979  G4int items_array_length=0;
980  char * buffer= new char[LINEBUFFSIZE];
981  size_t rflag = 0;
982 
983  while(items_array_length < SQ_Length)
984  {
985  rflag = std::fread(buffer, 2, 1, nested);
986  GetValue(buffer, item_GroupNumber);
987 
988  rflag = std::fread(buffer, 2, 1, nested);
989  GetValue(buffer, item_ElementNumber);
990 
991  rflag = std::fread(buffer, 4, 1, nested);
992  GetValue(buffer, item_Length);
993 
994  rflag = std::fread(buffer, item_Length, 1, nested);
995 
996  items_array_length= items_array_length+8+item_Length;
997  }
998 
999  delete [] buffer;
1000 
1001  if( SQ_Length>items_array_length )
1002  return 0;
1003  else
1004  return 1;
1005  if (rflag) return 1;
1006 }
1007 
1008 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1009 
1011 {
1012  // VARIABLES
1013  unsigned short item_GroupNumber;
1014  unsigned short item_ElementNumber;
1015  unsigned int item_Length;
1016  char * buffer= new char[LINEBUFFSIZE];
1017  size_t rflag = 0;
1018 
1019  do
1020  {
1021  rflag = std::fread(buffer, 2, 1, nested);
1022  GetValue(buffer, item_GroupNumber);
1023 
1024  rflag = std::fread(buffer, 2, 1, nested);
1025  GetValue(buffer, item_ElementNumber);
1026 
1027  rflag = std::fread(buffer, 4, 1, nested);
1028  GetValue(buffer, item_Length);
1029 
1030  if(item_Length!=0xffffffff)
1031  rflag = std::fread(buffer, item_Length, 1, nested);
1032  else
1033  read_undefined_item(nested);
1034 
1035 
1036  } while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD
1037  || item_Length!=0);
1038 
1039  delete [] buffer;
1040  if (rflag) return;
1041 }
1042 
1043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1044 
1046 {
1047  // VARIABLES
1048  unsigned short item_GroupNumber;
1049  unsigned short item_ElementNumber;
1050  G4int item_Length; size_t rflag = 0;
1051  char *buffer= new char[LINEBUFFSIZE];
1052 
1053  do
1054  {
1055  rflag = std::fread(buffer, 2, 1, nested);
1056  GetValue(buffer, item_GroupNumber);
1057 
1058  rflag = std::fread(buffer, 2, 1, nested);
1059  GetValue(buffer, item_ElementNumber);
1060 
1061  rflag = std::fread(buffer, 4, 1, nested);
1062  GetValue(buffer, item_Length);
1063 
1064 
1065  if(item_Length!=0)
1066  rflag = std::fread(buffer,item_Length,1,nested);
1067 
1068  }
1069  while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D
1070  || item_Length!=0);
1071 
1072  delete [] buffer;
1073  if (rflag) return;
1074 }
1075 
1076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1077 
1078 template <class Type>
1079 void DicomHandler::GetValue(char * _val, Type & _rval) {
1080 
1081 #if BYTE_ORDER == BIG_ENDIAN
1082  if(fLittleEndian) { // little endian
1083 #else // BYTE_ORDER == LITTLE_ENDIAN
1084  if(!fLittleEndian) { // big endian
1085 #endif
1086  const int SIZE = sizeof(_rval);
1087  char ctemp;
1088  for(int i = 0; i < SIZE/2; i++) {
1089  ctemp = _val[i];
1090  _val[i] = _val[SIZE - 1 - i];
1091  _val[SIZE - 1 - i] = ctemp;
1092  }
1093  }
1094  _rval = *(Type *)_val;
1095  }
1096 
1097  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1098 
G4double fSliceLocation
G4String driverFile
void GetValue(char *, Type &)
void SetMinZ(const G4double &val)
void SetMinX(const G4double &val)
G4int fMinPixelValue
G4String name
Definition: TRTMaterials.hh:40
bool readCalibration
void SetNoVoxelY(const G4int &val)
float G4float
Definition: G4Types.hh:77
void read_undefined_nested(FILE *)
void GetInformation(G4int &, char *)
void ReadCalibration()
void SetNoVoxelX(const G4int &val)
#define buffer
Definition: xmlparse.cc:611
short fBitAllocated
Definition of the DicomPhantomZSliceMerged class.
const int DATABUFFSIZE
Definition: DicomHandler.hh:97
void SetMinY(const G4double &val)
int G4int
Definition: G4Types.hh:78
DicomPhantomZSliceMerged * mergedSlices
G4int fRescaleSlope
static DicomHandler * fgInstance
DicomHandler.cc :
Definition: DicomHandler.hh:95
G4double density
Definition: TRTMaterials.hh:39
static DicomHandler * Instance()
Definition: DicomHandler.cc:72
G4bool fLittleEndian
G4int ** fTab
G4GLOB_DLL std::ostream G4cout
void SetMaxX(const G4double &val)
void SetSliceLocation(const G4double &val)
bool G4bool
Definition: G4Types.hh:79
DicomPhantomZSliceHeader class.
const int FILENAMESIZE
Definition: DicomHandler.hh:99
G4double fPixelSpacingX
unsigned int GetMaterialIndex(G4float density)
G4int ReadFile(FILE *, char *)
void StoreData(std::ofstream &foutG4DCM)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const char * data() const
G4bool fImplicitEndian
G4String ct2densityFile
G4double fSliceThickness
G4int ReadData(FILE *, char *)
G4double fPixelSpacingY
short fPixelRepresentation
Definition of the DicomHandler class.
G4double * valuedensity
G4float Pixel2density(G4int pixel)
#define G4endl
Definition: G4ios.hh:61
std::map< G4float, G4String > fMaterialIndices
void AddZSlice(DicomPhantomZSliceHeader *val)
void CheckFileFormat()
G4int fMaxPixelValue
Definition of the DicomPhantomZSliceHeader class.
void ReadMaterialIndices(std::ifstream &finData)
double G4double
Definition: G4Types.hh:76
void SetMaxZ(const G4double &val)
const int LINEBUFFSIZE
Definition: DicomHandler.hh:98
G4int read_defined_nested(FILE *, G4int)
void read_undefined_item(FILE *)
G4int fRescaleIntercept
void SetNoVoxelZ(const G4int &val)
G4double * valueCT
short fCompression
void AddMaterial(const G4String &val)
void SetMaxY(const G4double &val)