| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include "DicomFileMgr.hh"
27
28 #include "DicomFileCT.hh"
29 #include "DicomFilePET.hh"
30 #include "DicomFilePlan.hh"
31 #include "DicomFileStructure.hh"
32 #include "dcmtk/dcmdata/dcdeftag.h"
33
34 #include "G4UIcommand.hh"
35 #include "G4tgrFileIn.hh"
36
37 DicomFileMgr* DicomFileMgr::theInstance = 0;
38 int DicomFileMgr::verbose = 1;
39
40 //....oooOO0OOooo........oooOO0OOooo........oo
41 DicomFileMgr* DicomFileMgr::GetInstance()
42 {
43 if (!theInstance) {
44 theInstance = new DicomFileMgr;
45 }
46 return theInstance;
47 }
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50 DicomFileMgr::DicomFileMgr()
51 {
52 fCompression = 1.;
53 theCTFileAll = 0;
54 theStructureNCheck = 4;
55 theStructureNMaxROI = 100;
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59 void DicomFileMgr::Convert(G4String fileName)
60 {
61 G4tgrFileIn fin = G4tgrFileIn::GetInstance(f
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 }
72 else if (wl[0] == ":FILE") {
73 CheckNColumns(wl, 2);
74 G4cout << "@@@@@@@ Reading FILE: " << wl
75 AddFile(wl[1]);
76 }
77 else if (wl[0] == ":FILE_OUT") {
78 CheckNColumns(wl, 2);
79 theFileOutName = wl[1];
80 }
81 else if (wl[0] == ":MATE_DENS") {
82 CheckNColumns(wl, 3);
83 AddMaterialDensity(wl);
84 }
85 else if (wl[0] == ":MATE") {
86 CheckNColumns(wl, 3);
87 AddMaterial(wl);
88 }
89 else if (wl[0] == ":CT2D") {
90 CheckNColumns(wl, 3);
91 AddCT2Density(wl);
92 }
93 else {
94 G4Exception("DICOM2G4", "Wrong argument"
95 G4String("UNKNOWN TAG IN FIL
96 }
97 }
98
99 //@@@@@@ Process files
100 ProcessFiles();
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oo
104 void DicomFileMgr::CheckNColumns(std::vector<G
105 {
106 if (wl.size() != vsizeTh) {
107 G4cerr << " Reading line " << G4endl;
108 for (size_t ii = 0; ii < wl.size(); ii++)
109 G4cerr << wl[ii] << " ";
110 }
111 G4cerr << G4endl;
112 G4Exception("DICOM2G4", "D2G0010", FatalEr
113 ("Wrong number of columns in l
114 + std::to_string(vsizeTh))
115 .c_str());
116 }
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oo
120 void DicomFileMgr::SetCompression(G4String fCo
121 {
122 fCompression = G4UIcommand::ConvertToDouble(
123 }
124
125 //....oooOO0OOooo........oooOO0OOooo........oo
126 void DicomFileMgr::AddFile(G4String fileName)
127 {
128 DcmFileFormat dfile;
129 if (!(dfile.loadFile(fileName.c_str())).good
130 G4Exception("DicomHandler::ReadFile", "dfi
131 ("Error reading file " + fileN
132 }
133 DcmDataset* dset = dfile.getDataset();
134
135 OFString dModality;
136 if (!dset->findAndGetOFString(DCM_Modality,
137 G4Exception("DicomHandler::ReadData ", "",
138 }
139
140 if (dModality == "CT") {
141 DicomFileCT* df = new DicomFileCT(dset);
142 df->ReadData();
143 df->SetFileName(fileName);
144 // reorder by location
145 theCTFiles[df->GetMaxZ()] = df;
146 }
147 else if (dModality == "RTSTRUCT") {
148 DicomFileStructure* df = new DicomFileStru
149 df->ReadData();
150 df->SetFileName(fileName);
151 // theFiles.insert(msd::value_type(dMod
152 theStructFiles.push_back(df);
153 }
154 else if (dModality == "RTPLAN") {
155 DicomFilePlan* df = new DicomFilePlan(dset
156 df->ReadData();
157 df->SetFileName(fileName);
158 // theFiles.insert(msd::value_type(dMod
159 thePlanFiles.push_back(df);
160 }
161 else if (dModality == "PT") {
162 DicomFilePET* df = new DicomFilePET(dset);
163 df->ReadData();
164 df->SetFileName(fileName);
165 // theFiles.insert(msd::value_type(dMod
166 thePETFiles[df->GetMaxZ()] = df;
167 // thePETFiles.push_back(df);
168 }
169 else {
170 G4Exception(
171 "DicomFileMgr::AddFIle()", "DFM001", Fat
172 (G4String("File is not of type CT or RTS
173 }
174 }
175
176 //....oooOO0OOooo........oooOO0OOooo........oo
177 void DicomFileMgr::AddMaterial(std::vector<G4S
178 {
179 if (theMaterials.size() > 0 && bMaterialsDen
180 G4Exception(
181 "DicomFileMgr::AddMaterial", "DFM005", F
182 "Trying to add a Material with :MATE and
183 }
184 bMaterialsDensity = false;
185 theMaterials[G4UIcommand::ConvertToDouble(wl
186 }
187
188 //....oooOO0OOooo........oooOO0OOooo........oo
189 void DicomFileMgr::AddMaterialDensity(std::vec
190 {
191 if (theMaterialsDensity.size() > 0 && !bMate
192 G4Exception(
193 "DicomFileMgr::AddMaterial", "DFM005", F
194 "Trying to add a Material with :MATE and
195 }
196 bMaterialsDensity = true;
197 theMaterialsDensity[G4UIcommand::ConvertToDo
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oo
201 void DicomFileMgr::AddCT2Density(std::vector<G
202 {
203 theCT2Density[G4UIcommand::ConvertToInt(wl[1
204 G4cout << this << " AddCT2density " << theCT
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oo
208 G4double DicomFileMgr::Hounsfield2density(Uint
209 {
210 if (theCT2Density.size() == 0) {
211 G4Exception("Hounsfield2density", "DCM004"
212 }
213 std::map<G4int, G4double>::const_iterator it
214 G4int minHval = (*ite).first;
215 if (G4int(Hval) < minHval) {
216 G4Exception("DicomHandler::Hounsfield2dens
217 ("Hval value too small, change
218 + std::to_string(minHval))
219 .c_str());
220 }
221
222 ite = theCT2Density.end();
223 ite--;
224 G4int maxHval = (*ite).first;
225 if (G4int(Hval) > maxHval) {
226 G4Exception("DicomHandler::Hval2density",
227 ("Hval value too big, change C
228 + std::to_string(maxHval))
229 .c_str());
230 }
231
232 G4float density = -1.;
233 G4double deltaCT = 0;
234 G4double deltaDensity = 0;
235
236 ite = theCT2Density.upper_bound(Hval);
237 std::map<G4int, G4double>::const_iterator it
238 itePrev--;
239
240 deltaCT = (*ite).first - (*itePrev).first;
241 deltaDensity = (*ite).second - (*itePrev).se
242
243 // interpolating linearly
244 density = (*ite).second - (((*ite).first - H
245
246 if (density < 0.) {
247 G4Exception("DicomFileMgr::Hounsfiled2Dens
248 G4String("@@@ Error negative d
249 + " from HV = " + std
250 .c_str());
251 }
252
253 // G4cout << " Hval2density " << Hval << "
254 return density;
255 }
256
257 //....oooOO0OOooo........oooOO0OOooo........oo
258 size_t DicomFileMgr::GetMaterialIndex(G4double
259 {
260 std::map<G4double, G4String>::iterator ite =
261 if (ite == theMaterials.end()) {
262 ite--;
263 G4Exception("DicomFileMgr::GetMaterialInde
264 ("Hounsfiled value too big, ch
265 + std::to_string((*ite).first
266 .c_str());
267 }
268
269 size_t dist = std::distance(theMaterials.beg
270
271 return dist;
272 }
273
274 //....oooOO0OOooo........oooOO0OOooo........oo
275 size_t DicomFileMgr::GetMaterialIndexByDensity
276 {
277 std::map<G4double, G4String>::iterator ite =
278 if (ite == theMaterialsDensity.end()) {
279 ite--;
280 G4Exception("DicomFileMgr::GetMaterialInde
281 ("Density too big, change inpu
282 + std::to_string((*ite).first
283 .c_str());
284 }
285
286 size_t dist = std::distance(theMaterialsDens
287
288 return dist;
289 }
290
291 //....oooOO0OOooo........oooOO0OOooo........oo
292 void DicomFileMgr::ProcessFiles()
293 {
294 if (theCTFiles.size() == 0) {
295 G4Exception("CheckCTSlices", "DCM004", Jus
296 }
297 else {
298 CheckCTSlices();
299
300 BuildCTMaterials();
301
302 MergeCTFiles();
303 }
304
305 G4cout << " PROCESSING PET FILES " << thePET
306 if (thePETFiles.size() != 0) {
307 CheckPETSlices();
308
309 BuildPETActivities();
310
311 MergePETFiles();
312 }
313
314 DumpToTextFile();
315 }
316
317 //....oooOO0OOooo........oooOO0OOooo........oo
318 void DicomFileMgr::CheckCTSlices()
319 {
320 size_t nSlices = theCTFiles.size();
321 G4cout << " DicomFileMgr::Checking CT slices
322
323 G4bool uniformSliceThickness = true;
324
325 if (nSlices > 1) {
326 if (nSlices == 2) {
327 mdct::const_iterator ite = theCTFiles.be
328 DicomFileCT* one = (*ite).second;
329 ite++;
330 DicomFileCT* two = (*ite).second;
331
332 G4double real_distance = (two->GetLocati
333
334 if (one->GetMaxZ() != two->GetMinZ()) {
335 one->SetMaxZ(one->GetLocation() + real
336 two->SetMinZ(two->GetLocation() - real
337 // one->SetMinZ(one->GetLocation()-rea
338 // two->SetMaxZ(two->GetLocation()+rea
339 if (uniformSliceThickness) {
340 one->SetMinZ(one->GetLocation() - re
341 two->SetMaxZ(two->GetLocation() + re
342 }
343 }
344 }
345 else {
346 mdct::iterator ite0 = theCTFiles.begin()
347 mdct::iterator ite1 = ite0;
348 ite1++;
349 mdct::iterator ite2 = ite1;
350 ite2++;
351 for (; ite2 != theCTFiles.end(); ++ite0,
352 DicomFileCT* prev = (DicomFileCT*)((*i
353 DicomFileCT* slice = (DicomFileCT*)((*
354 DicomFileCT* next = (DicomFileCT*)((*i
355 G4double real_up_distance = (next->Get
356 G4double real_down_distance = (slice->
357 G4double real_distance = real_up_dista
358 G4double stated_distance = slice->GetM
359
360 if (std::fabs(real_distance - stated_d
361 unsigned int sliceNum = std::distanc
362 G4cerr << "\tDicomFileMgr::CheckCTSl
363 << "]: Distance between this
364 << " <> Slice width = " << st
365 << prev->GetLocation() << " :
366 << next->GetLocation() << " D
367 << G4endl;
368 G4cerr << "!! WARNING: Geant4 will r
369 << "lower and upper slice " <
370
371 slice->SetMinZ(slice->GetLocation()
372 slice->SetMaxZ(slice->GetLocation()
373
374 if (ite0 == theCTFiles.begin()) {
375 prev->SetMaxZ(slice->GetMinZ());
376 // Using below would make all slic
377 // prev->SetMinZ(prev->GetLocation
378 if (uniformSliceThickness) {
379 prev->SetMinZ(prev->GetLocation(
380 }
381 }
382 if (static_cast<unsigned int>(std::d
383 next->SetMinZ(slice->GetMaxZ());
384 // Using below would make all slic
385 // next->SetMaxZ(next->GetLocation
386 if (uniformSliceThickness) {
387 next->SetMaxZ(next->GetLocation(
388 }
389 }
390 }
391 }
392 }
393 }
394 }
395
396 //....oooOO0OOooo........oooOO0OOooo........oo
397 void DicomFileMgr::CheckPETSlices()
398 {
399 size_t nSlices = thePETFiles.size();
400 G4cout << " DicomFileMgr::Checking PET slice
401
402 G4bool uniformSliceThickness = true;
403
404 if (nSlices > 1) {
405 if (nSlices == 2) {
406 mdpet::const_iterator ite = thePETFiles.
407 DicomFilePET* one = (*ite).second;
408 ite++;
409 DicomFilePET* two = (*ite).second;
410
411 G4double real_distance = (two->GetLocati
412
413 if (one->GetMaxZ() != two->GetMinZ()) {
414 one->SetMaxZ(one->GetLocation() + real
415 two->SetMinZ(two->GetLocation() - real
416 // one->SetMinZ(one->GetLocation()-rea
417 // two->SetMaxZ(two->GetLocation()+rea
418 if (uniformSliceThickness) {
419 one->SetMinZ(one->GetLocation() - re
420 two->SetMaxZ(two->GetLocation() + re
421 }
422 }
423 }
424 else {
425 mdpet::iterator ite0 = thePETFiles.begin
426 mdpet::iterator ite1 = ite0;
427 ite1++;
428 mdpet::iterator ite2 = ite1;
429 ite2++;
430 for (; ite2 != thePETFiles.end(); ++ite0
431 DicomFilePET* prev = (DicomFilePET*)((
432 DicomFilePET* slice = (DicomFilePET*)(
433 DicomFilePET* next = (DicomFilePET*)((
434 G4double real_up_distance = (next->Get
435 G4double real_down_distance = (slice->
436 G4double real_distance = real_up_dista
437 G4double stated_distance = slice->GetM
438
439 if (std::fabs(real_distance - stated_d
440 unsigned int sliceNum = std::distanc
441 G4cerr << "\tDicomFileMgr::CheckPETS
442 << "]: Distance between this
443 << " <> Slice width = " << st
444 << prev->GetLocation() << " :
445 << next->GetLocation() << " D
446 << G4endl;
447 G4cerr << "!! WARNING: Geant4 will r
448 << "lower and upper slice " <
449
450 slice->SetMinZ(slice->GetLocation()
451 slice->SetMaxZ(slice->GetLocation()
452
453 if (ite0 == thePETFiles.begin()) {
454 prev->SetMaxZ(slice->GetMinZ());
455 // Using below would make all slic
456 // prev->SetMinZ(prev->GetLocation
457 if (uniformSliceThickness) {
458 prev->SetMinZ(prev->GetLocation(
459 }
460 }
461 if (static_cast<unsigned int>(std::d
462 next->SetMinZ(slice->GetMaxZ());
463 // Using below would make all slic
464 // next->SetMaxZ(next->GetLocation
465 if (uniformSliceThickness) {
466 next->SetMaxZ(next->GetLocation(
467 }
468 }
469 }
470 }
471 }
472 }
473 }
474
475 //....oooOO0OOooo........oooOO0OOooo........oo
476 void DicomFileMgr::BuildCTMaterials()
477 {
478 G4cout << " DicomFileMgr::Building Materials
479 mdct::const_iterator ite = theCTFiles.begin(
480 for (; ite != theCTFiles.end(); ite++) {
481 (*ite).second->BuildMaterials();
482 }
483 }
484
485 //....oooOO0OOooo........oooOO0OOooo........oo
486 void DicomFileMgr::BuildPETActivities()
487 {
488 G4cout << " DicomFileMgr::Building PETData "
489 mdpet::const_iterator ite = thePETFiles.begi
490 for (; ite != thePETFiles.end(); ite++) {
491 (*ite).second->BuildActivities();
492 }
493 }
494
495 //....oooOO0OOooo........oooOO0OOooo........oo
496 void DicomFileMgr::MergeCTFiles()
497 {
498 G4cout << " DicomFileMgr::Merging CT Files "
499 mdct::const_iterator ite = theCTFiles.begin(
500 theCTFileAll = new DicomFileCT(*((*ite).seco
501 ite++;
502 for (; ite != theCTFiles.end(); ite++) {
503 (*theCTFileAll) += *((*ite).second);
504 }
505 }
506
507 //....oooOO0OOooo........oooOO0OOooo........oo
508 void DicomFileMgr::MergePETFiles()
509 {
510 G4cout << " DicomFileMgr::Merging PET Files
511 mdpet::const_iterator ite = thePETFiles.begi
512 thePETFileAll = new DicomFilePET(*((*ite).se
513 ite++;
514 for (; ite != thePETFiles.end(); ite++) {
515 (*thePETFileAll) += *((*ite).second);
516 }
517 }
518
519 //....oooOO0OOooo........oooOO0OOooo........oo
520 void DicomFileMgr::DumpToTextFile()
521 {
522 G4cout << " DicomFileMgr::Dumping To Text Fi
523 if (theCTFiles.size() != 0) {
524 std::ofstream fout(theFileOutName);
525
526 if (!bMaterialsDensity) {
527 fout << theMaterials.size() << std::endl
528 std::map<G4double, G4String>::const_iter
529 G4int ii = 0;
530 for (ite = theMaterials.begin(); ite !=
531 fout << ii << " \"" << (*ite).second <
532 }
533 }
534 else {
535 fout << theMaterialsDensity.size() << st
536 std::map<G4double, G4String>::const_iter
537 G4int ii = 0;
538 for (ite = theMaterialsDensity.begin();
539 fout << ii << " \"" << (*ite).second <
540 }
541 }
542
543 theCTFileAll->DumpHeaderToTextFile(fout);
544 for (mdct::const_iterator itect = theCTFil
545 (*itect).second->DumpMateIDsToTextFile(f
546 }
547 for (mdct::const_iterator itect = theCTFil
548 (*itect).second->DumpDensitiesToTextFile
549 }
550 for (mdct::const_iterator itect = theCTFil
551 (*itect).second->BuildStructureIDs();
552 (*itect).second->DumpStructureIDsToTextF
553 }
554
555 std::vector<DicomFileStructure*> dfs = Get
556 for (size_t i1 = 0; i1 < dfs.size(); i1++)
557 std::vector<DicomROI*> rois = dfs[i1]->G
558 for (size_t i2 = 0; i2 < rois.size(); i2
559 fout << rois[i2]->GetNumber() + 1 << "
560 }
561 }
562 }
563
564 if (thePETFiles.size() != 0) {
565 std::ofstream fout(theFileOutName);
566
567 thePETFileAll->DumpHeaderToTextFile(fout);
568 for (mdpet::const_iterator itect = thePETF
569 (*itect).second->DumpActivitiesToTextFil
570 }
571 }
572
573 for (size_t i1 = 0; i1 < thePlanFiles.size()
574 thePlanFiles[i1]->DumpToFile();
575 }
576 }
577