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