ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomFileMgr.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DicomFileMgr.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 #include "DicomFileMgr.hh"
27 
28 #include "DicomFileCT.hh"
29 #include "DicomFilePET.hh"
30 #include "DicomFileStructure.hh"
31 #include "DicomFilePlan.hh"
32 
33 #include "dcmtk/dcmdata/dcdeftag.h"
34 #include "G4tgrFileIn.hh"
35 #include "G4UIcommand.hh"
36 
38 int DicomFileMgr::verbose = 1;
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 {
43  if( !theInstance ) {
45  }
46  return theInstance;
47 }
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 {
52  fCompression = 1.;
53  theCTFileAll = 0;
55  theStructureNMaxROI = 100;
56 }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 {
62  std::vector<G4String> wl;
63  // Read each file in file list
64  theFileOutName = "test.g4dcm";
65  int ii;
66  for( ii = 0;; ii++) {
67  if( ! fin.GetWordsInLine(wl) ) break;
68  if( wl[0] == ":COMPRESSION" ) {
69  CheckNColumns(wl,2);
70  SetCompression(wl[1]);
71  } else if( wl[0] == ":FILE" ) {
72  CheckNColumns(wl,2);
73  G4cout << "@@@@@@@ Reading FILE: " << wl[1] << G4endl;
74  AddFile(wl[1]);
75  } else if( wl[0] == ":FILE_OUT" ) {
76  CheckNColumns(wl,2);
77  theFileOutName = wl[1];
78  } else if( wl[0] == ":MATE_DENS" ) {
79  CheckNColumns(wl,3);
81  } else if( wl[0] == ":MATE" ) {
82  CheckNColumns(wl,3);
83  AddMaterial(wl);
84  } else if( wl[0] == ":CT2D" ) {
85  CheckNColumns(wl,3);
86  AddCT2Density(wl);
87  } else {
88  G4Exception("DICOM2G4",
89  "Wrong argument",
91  G4String("UNKNOWN TAG IN FILE "+wl[0]).c_str());
92  }
93 
94  }
95 
96 
97  //@@@@@@ Process files
98  ProcessFiles();
99 
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 void DicomFileMgr::CheckNColumns(std::vector<G4String> wl, size_t vsizeTh )
104 {
105  if( wl.size() != vsizeTh ) {
106  G4cerr << " Reading line " << G4endl;
107  for( size_t ii = 0; ii < wl.size(); ii++){
108  G4cerr << wl[ii] << " ";
109  }
110  G4cerr << G4endl;
111  G4Exception("DICOM2G4",
112  "D2G0010",
114  ("Wrong number of columns in line " + std::to_string(wl.size()) + " <> "
115  + std::to_string(vsizeTh)).c_str());
116  }
117 
118 }
119 
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123 {
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 {
130  DcmFileFormat dfile;
131  if( ! (dfile.loadFile(fileName.c_str())).good() ) {
132  G4Exception("DicomHandler::ReadFile",
133  "dfile.loadFile",
135  ("Error reading file " + fileName).c_str());
136  }
137  DcmDataset* dset = dfile.getDataset();
138 
139  OFString dModality;
140  if( !dset->findAndGetOFString(DCM_Modality,dModality).good() ) {
141  G4Exception("DicomHandler::ReadData ",
142  "",
144  " Have not read Modality");
145  }
146 
147  if( dModality == "CT" ) {
148  DicomFileCT* df = new DicomFileCT(dset);
149  df->ReadData();
150  df->SetFileName( fileName );
151  // reorder by location
152  theCTFiles[df->GetMaxZ()] = df;
153  } else if( dModality == "RTSTRUCT" ) {
154  DicomFileStructure* df = new DicomFileStructure(dset);
155  df->ReadData();
156  df->SetFileName( fileName );
157  // theFiles.insert(msd::value_type(dModality,df));
158  theStructFiles.push_back(df);
159  } else if( dModality == "RTPLAN" ) {
160  DicomFilePlan* df = new DicomFilePlan(dset);
161  df->ReadData();
162  df->SetFileName( fileName );
163  // theFiles.insert(msd::value_type(dModality,df));
164  thePlanFiles.push_back(df);
165  } else if( dModality == "PT" ) {
166  DicomFilePET* df = new DicomFilePET(dset);
167  df->ReadData();
168  df->SetFileName( fileName );
169  // theFiles.insert(msd::value_type(dModality,df));
170  thePETFiles[df->GetMaxZ()] = df;
171  // thePETFiles.push_back(df);
172  } else {
173  G4Exception("DicomFileMgr::AddFIle()",
174  "DFM001",
176  (G4String("File is not of type CT or RTSTRUCT or RTPLAN, but: ")
177  + dModality).c_str());
178  }
179 
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183 void DicomFileMgr::AddMaterial( std::vector<G4String> wl )
184 {
185  if( theMaterials.size() > 0 && bMaterialsDensity ) {
186  G4Exception("DicomFileMgr::AddMaterial",
187  "DFM005",
189  "Trying to add a Material with :MATE and another with :MATE_DENS, check your input file");
190  }
191  bMaterialsDensity = false;
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 void DicomFileMgr::AddMaterialDensity( std::vector<G4String> wl )
197 {
198  if( theMaterialsDensity.size() > 0 && !bMaterialsDensity ) {
199  G4Exception("DicomFileMgr::AddMaterial",
200  "DFM005",
202  "Trying to add a Material with :MATE and another with :MATE_DENS, check your input file");
203  }
204  bMaterialsDensity = true;
206 }
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 void DicomFileMgr::AddCT2Density( std::vector<G4String> wl)
210 {
212  G4cout << this << " AddCT2density " << theCT2Density.size() << G4endl;//GDEB
213 
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 {
219  if( theCT2Density.size() == 0 ) {
220  G4Exception("Hounsfield2density",
221  "DCM004",
223  "No :CT2D line in input file");
224  }
225  std::map<G4int,G4double>::const_iterator ite = theCT2Density.begin();
226  G4int minHval = (*ite).first;
227  if( G4int(Hval) < minHval ) {
228  G4Exception("DicomHandler::Hounsfield2density",
229  "",
231  ("Hval value too small, change input file "+std::to_string(Hval) + " < "
232  + std::to_string(minHval)).c_str());
233  }
234 
235  ite = theCT2Density.end(); ite--;
236  G4int maxHval = (*ite).first;
237  if( G4int(Hval) > maxHval ) {
238  G4Exception("DicomHandler::Hval2density",
239  "",
241  ("Hval value too big, change CT2Density.dat file "+std::to_string(Hval) + " > "
242  + std::to_string(maxHval)).c_str());
243  }
244 
245  G4float density = -1.;
246  G4double deltaCT = 0;
247  G4double deltaDensity = 0;
248 
249  ite = theCT2Density.upper_bound(Hval);
250  std::map<G4int,G4double>::const_iterator itePrev = ite; itePrev--;
251 
252  deltaCT = (*ite).first - (*itePrev).first;
253  deltaDensity = (*ite).second - (*itePrev).second;
254 
255  // interpolating linearly
256  density = (*ite).second - (((*ite).first-Hval)*deltaDensity/deltaCT );
257 
258  if(density < 0.) {
259  G4Exception("DicomFileMgr::Hounsfiled2Density",
260  "DFM002",
262  G4String("@@@ Error negative density = " + std::to_string(density) + " from HV = "
263  + std::to_string(Hval)).c_str());
264  }
265 
266  // G4cout << " Hval2density " << Hval << " -> " << density << G4endl;//GDEB
267  return density;
268 
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273 {
274  std::map<G4double,G4String>::iterator ite = theMaterials.upper_bound(Hval);
275  if( ite == theMaterials.end() ) {
276  ite--;
277  G4Exception("DicomFileMgr::GetMaterialIndex",
278  "DFM004",
280  ("Hounsfiled value too big, change input file "+std::to_string(Hval) + " > "
281  + std::to_string((*ite).first)).c_str());
282  }
283 
284  size_t dist = std::distance( theMaterials.begin(), ite );
285 
286  return dist;
287 
288 }
289 
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
293 {
294  std::map<G4double,G4String>::iterator ite = theMaterialsDensity.upper_bound(density);
295  if( ite == theMaterialsDensity.end() ) {
296  ite--;
297  G4Exception("DicomFileMgr::GetMaterialIndexByDensity",
298  "DFM003",
300  ("Density too big, change input file "+std::to_string(density) + " > "
301  + std::to_string((*ite).first)).c_str());
302  }
303 
304  size_t dist = std::distance( theMaterialsDensity.begin(), ite );
305 
306  return dist;
307 
308 }
309 
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 {
313  if( theCTFiles.size() == 0 ) {
314  G4Exception("CheckCTSlices",
315  "DCM004",
316  JustWarning,
317  "No :FILE of type CT in input file");
318  } else {
319 
320  CheckCTSlices();
321 
323 
324  MergeCTFiles();
325 
326  }
327 
328  G4cout << " PROCESSING PET FILES " << thePETFiles.size() << G4endl; //GDEB
329  if( thePETFiles.size() != 0 ) {
330 
331  CheckPETSlices();
332 
334 
335  MergePETFiles();
336 
337  }
338 
339  DumpToTextFile();
340 
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345 {
346  size_t nSlices = theCTFiles.size();
347  G4cout << " DicomFileMgr::Checking CT slices: " << nSlices << G4endl;
348 
349  G4bool uniformSliceThickness = true;
350 
351  if(nSlices > 1) {
352  if(nSlices == 2) {
353  mdct::const_iterator ite = theCTFiles.begin();
354  DicomFileCT* one = (*ite).second;
355  ite++;
356  DicomFileCT* two = (*ite).second;
357 
358  G4double real_distance = (two->GetLocation()-one->GetLocation())/2.;
359 
360  if(one->GetMaxZ() != two->GetMinZ()) {
361  one->SetMaxZ(one->GetLocation()+real_distance);
362  two->SetMinZ(two->GetLocation()-real_distance);
363  //one->SetMinZ(one->GetLocation()-real_distance);
364  //two->SetMaxZ(two->GetLocation()+real_distance);
365  if(uniformSliceThickness) {
366  one->SetMinZ(one->GetLocation()-real_distance);
367  two->SetMaxZ(two->GetLocation()+real_distance);
368  }
369  }
370  } else {
371  mdct::iterator ite0 = theCTFiles.begin();
372  mdct::iterator ite1 = ite0; ite1++;
373  mdct::iterator ite2 = ite1; ite2++;
374  for(; ite2 != theCTFiles.end(); ++ite0, ++ite1, ++ite2) {
375  DicomFileCT* prev = (DicomFileCT*)((*ite0).second);
376  DicomFileCT* slice = (DicomFileCT*)((*ite1).second);
377  DicomFileCT* next = (DicomFileCT*)((*ite2).second);
378  G4double real_up_distance = (next->GetLocation() -
379  slice->GetLocation())/2.;
380  G4double real_down_distance = (slice->GetLocation() -
381  prev->GetLocation())/2.;
382  G4double real_distance = real_up_distance + real_down_distance;
383  G4double stated_distance = slice->GetMaxZ()-slice->GetMinZ();
384 
385  if(std::fabs(real_distance - stated_distance) > 1.E-9) {
386  unsigned int sliceNum = std::distance(theCTFiles.begin(),ite1);
387  G4cerr << "\tDicomFileMgr::CheckCTSlices - Slice Distance Error in slice [" << sliceNum
388  << "]: Distance between this slice and slices up and down = "
389  << real_distance
390  << " <> Slice width = " << stated_distance
391  << " Slice locations " <<prev->GetLocation() << " : " << slice->GetLocation()
392  << " : " << next->GetLocation()
393  << " DIFFERENCE= " << real_distance - stated_distance
394  << G4endl;
395  G4cerr << "!! WARNING: Geant4 will reset slice width so that it extends between "
396  << "lower and upper slice " << G4endl;
397 
398  slice->SetMinZ(slice->GetLocation()-real_down_distance);
399  slice->SetMaxZ(slice->GetLocation()+real_up_distance);
400 
401  if(ite0 == theCTFiles.begin()) {
402  prev->SetMaxZ(slice->GetMinZ());
403  // Using below would make all slice same thickness
404  //prev->SetMinZ(prev->GetLocation()-real_min_distance);
405  if(uniformSliceThickness) {
406  prev->SetMinZ(prev->GetLocation()-real_down_distance);
407  }
408  }
409  if(static_cast<unsigned int>(std::distance(theCTFiles.begin(),ite2)+1)==
410  nSlices) {
411  next->SetMinZ(slice->GetMaxZ());
412  // Using below would make all slice same thickness
413  //next->SetMaxZ(next->GetLocation()+real_max_distance);
414  if(uniformSliceThickness) {
415  next->SetMaxZ(next->GetLocation()+real_up_distance); }
416  }
417  }
418  }
419  }
420  }
421 
422 }
423 
424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
426 {
427  size_t nSlices = thePETFiles.size();
428  G4cout << " DicomFileMgr::Checking PET slices: " << nSlices << G4endl;
429 
430  G4bool uniformSliceThickness = true;
431 
432  if(nSlices > 1) {
433  if(nSlices == 2) {
434  mdpet::const_iterator ite = thePETFiles.begin();
435  DicomFilePET* one = (*ite).second;
436  ite++;
437  DicomFilePET* two = (*ite).second;
438 
439  G4double real_distance = (two->GetLocation()-one->GetLocation())/2.;
440 
441  if(one->GetMaxZ() != two->GetMinZ()) {
442  one->SetMaxZ(one->GetLocation()+real_distance);
443  two->SetMinZ(two->GetLocation()-real_distance);
444  //one->SetMinZ(one->GetLocation()-real_distance);
445  //two->SetMaxZ(two->GetLocation()+real_distance);
446  if(uniformSliceThickness) {
447  one->SetMinZ(one->GetLocation()-real_distance);
448  two->SetMaxZ(two->GetLocation()+real_distance);
449  }
450  }
451  } else {
452  mdpet::iterator ite0 = thePETFiles.begin();
453  mdpet::iterator ite1 = ite0; ite1++;
454  mdpet::iterator ite2 = ite1; ite2++;
455  for(; ite2 != thePETFiles.end(); ++ite0, ++ite1, ++ite2) {
456  DicomFilePET* prev = (DicomFilePET*)((*ite0).second);
457  DicomFilePET* slice = (DicomFilePET*)((*ite1).second);
458  DicomFilePET* next = (DicomFilePET*)((*ite2).second);
459  G4double real_up_distance = (next->GetLocation() -
460  slice->GetLocation())/2.;
461  G4double real_down_distance = (slice->GetLocation() -
462  prev->GetLocation())/2.;
463  G4double real_distance = real_up_distance + real_down_distance;
464  G4double stated_distance = slice->GetMaxZ()-slice->GetMinZ();
465 
466  if(std::fabs(real_distance - stated_distance) > 1.E-9) {
467  unsigned int sliceNum = std::distance(thePETFiles.begin(),ite1);
468  G4cerr << "\tDicomFileMgr::CheckPETSlices - Slice Distance Error in slice [" << sliceNum
469  << "]: Distance between this slice and slices up and down = "
470  << real_distance
471  << " <> Slice width = " << stated_distance
472  << " Slice locations " <<prev->GetLocation() << " : " << slice->GetLocation()
473  << " : " << next->GetLocation()
474  << " DIFFERENCE= " << real_distance - stated_distance
475  << G4endl;
476  G4cerr << "!! WARNING: Geant4 will reset slice width so that it extends between "
477  << "lower and upper slice " << G4endl;
478 
479  slice->SetMinZ(slice->GetLocation()-real_down_distance);
480  slice->SetMaxZ(slice->GetLocation()+real_up_distance);
481 
482  if(ite0 == thePETFiles.begin()) {
483  prev->SetMaxZ(slice->GetMinZ());
484  // Using below would make all slice same thickness
485  //prev->SetMinZ(prev->GetLocation()-real_min_distance);
486  if(uniformSliceThickness) {
487  prev->SetMinZ(prev->GetLocation()-real_down_distance);
488  }
489  }
490  if(static_cast<unsigned int>(std::distance(thePETFiles.begin(),ite2)+1)==
491  nSlices) {
492  next->SetMinZ(slice->GetMaxZ());
493  // Using below would make all slice same thickness
494  //next->SetMaxZ(next->GetLocation()+real_max_distance);
495  if(uniformSliceThickness) {
496  next->SetMaxZ(next->GetLocation()+real_up_distance); }
497  }
498  }
499  }
500  }
501  }
502 
503 }
504 
505 
506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
508 {
509  G4cout << " DicomFileMgr::Building Materials " << theCTFiles.size() << G4endl;//GDEB
510  mdct::const_iterator ite = theCTFiles.begin();
511  for( ; ite != theCTFiles.end(); ite++ ) {
512  (*ite).second->BuildMaterials();
513  }
514 }
515 
516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
518 {
519  G4cout << " DicomFileMgr::Building PETData " << thePETFiles.size() << G4endl;//GDEB
520  mdpet::const_iterator ite = thePETFiles.begin();
521  for( ; ite != thePETFiles.end(); ite++ ) {
522  (*ite).second->BuildActivities();
523  }
524 }
525 
526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528 {
529  G4cout << " DicomFileMgr::Merging CT Files " << theCTFiles.size() << G4endl;//GDEB
530  mdct::const_iterator ite = theCTFiles.begin();
531  theCTFileAll = new DicomFileCT( *((*ite).second) );
532  ite++;
533  for( ; ite != theCTFiles.end(); ite++ ) {
534  (*theCTFileAll) += *((*ite).second);
535  }
536 }
537 
538 
539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
541 {
542  G4cout << " DicomFileMgr::Merging PET Files " << thePETFiles.size() << G4endl;//GDEB
543  mdpet::const_iterator ite = thePETFiles.begin();
544  thePETFileAll = new DicomFilePET( *((*ite).second) );
545  ite++;
546  for( ; ite != thePETFiles.end(); ite++ ) {
547  (*thePETFileAll) += *((*ite).second);
548  }
549 }
550 
551 
552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
554 {
555  G4cout << " DicomFileMgr::Dumping To Text File " << G4endl; //GDEB
556  if( theCTFiles.size() != 0 ) {
557  std::ofstream fout(theFileOutName);
558 
559  if( !bMaterialsDensity ) {
560  fout << theMaterials.size() << std::endl;
561  std::map<G4double,G4String>::const_iterator ite;
562  G4int ii = 0;
563  for(ite = theMaterials.begin(); ite != theMaterials.end(); ite++, ii++){
564  fout << ii << " \"" << (*ite).second << "\"" << std::endl;
565  }
566  } else {
567  fout << theMaterialsDensity.size() << std::endl;
568  std::map<G4double,G4String>::const_iterator ite;
569  G4int ii = 0;
570  for(ite = theMaterialsDensity.begin(); ite != theMaterialsDensity.end(); ite++, ii++){
571  fout << ii << " \"" << (*ite).second << "\"" << std::endl;
572  }
573  }
574 
576  for( mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++ ) {
577  (*itect).second->DumpMateIDsToTextFile(fout);
578  }
579  for( mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++ ) {
580  (*itect).second->DumpDensitiesToTextFile(fout);
581  }
582  for( mdct::const_iterator itect = theCTFiles.begin(); itect != theCTFiles.end(); itect++ ) {
583  (*itect).second->BuildStructureIDs();
584  (*itect).second->DumpStructureIDsToTextFile(fout);
585  }
586 
587  std::vector<DicomFileStructure*> dfs = GetStructFiles();
588  for( size_t i1 = 0; i1 < dfs.size(); i1++ ){
589  std::vector<DicomROI*> rois = dfs[i1]->GetROIs();
590  for( size_t i2 = 0; i2 < rois.size(); i2++ ){
591  fout << rois[i2]->GetNumber()+1 << " \"" << rois[i2]->GetName() << "\"" <<G4endl;
592  }
593  }
594  }
595 
596  if( thePETFiles.size() != 0 ) {
597  std::ofstream fout(theFileOutName);
598 
600  for( mdpet::const_iterator itect = thePETFiles.begin(); itect != thePETFiles.end(); itect++ ) {
601  (*itect).second->DumpActivitiesToTextFile(fout);
602  }
603  }
604 
605  for( size_t i1 = 0; i1 < thePlanFiles.size(); i1++ ){
606  thePlanFiles[i1]->DumpToFile();
607  }
608 
609 }