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