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