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