| 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 "DicomFileCT.hh"
27
28 #include "DicomFileStructure.hh"
29 #include "DicomROI.hh"
30 #include "dcmtk/dcmdata/dcdeftag.h"
31 #include "dcmtk/dcmdata/dcfilefo.h"
32 #include "dcmtk/dcmdata/dcpixel.h"
33 #include "dcmtk/dcmdata/dcpixseq.h"
34 #include "dcmtk/dcmdata/dcpxitem.h"
35 #include "dcmtk/dcmrt/drtimage.h"
36
37 #include "G4GeometryTolerance.hh"
38
39 #include <set>
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42 DicomFileCT::DicomFileCT() {}
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45 DicomFileCT::DicomFileCT(DcmDataset* dset) : D
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48 void DicomFileCT::BuildMaterials()
49 {
50 G4int fCompress = theFileMgr->GetCompression
51 if (fNoVoxelsX % fCompress != 0 || fNoVoxels
52 G4Exception("DicompFileMgr.:BuildMaterials
53 ("Compression factor = " + std
54 + " has to be a divisor of Nu
55 + " and Y " + std::to_string(
56 .c_str());
57 }
58
59 // if( DicomVerb(debugVerb) ) G4cout << " B
60 double meanHV = 0.;
61 for (int ir = 0; ir < fNoVoxelsY; ir += fCom
62 for (int ic = 0; ic < fNoVoxelsX; ic += fC
63 meanHV = 0.;
64 int isumrMax = std::min(ir + fCompress,
65 int isumcMax = std::min(ic + fCompress,
66 for (int isumr = ir; isumr < isumrMax; i
67 for (int isumc = ic; isumc < isumcMax;
68 meanHV += fHounsfieldV[isumc + isumr
69 // G4cout << isumr << " " << isumc <
70 }
71 }
72 meanHV /= (isumrMax - ir) * (isumcMax -
73 G4double meanDens = theFileMgr->Hounsfie
74 // G4cout << ir << " " << ic << " F
75
76 fDensities.push_back(meanDens);
77 size_t mateID;
78 if (theFileMgr->IsMaterialsDensity()) {
79 mateID = theFileMgr->GetMaterialIndexB
80 }
81 else {
82 mateID = theFileMgr->GetMaterialIndex(
83 }
84 fMateIDs.push_back(mateID);
85 }
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90 void DicomFileCT::DumpMateIDsToTextFile(std::o
91 {
92 G4int fCompress = theFileMgr->GetCompression
93 if (DicomFileMgr::verbose >= warningVerb)
94 G4cout << fLocation << " DumpMateIDsToText
95 << G4endl;
96 for (int ir = 0; ir < fNoVoxelsY / fCompress
97 for (int ic = 0; ic < fNoVoxelsX / fCompre
98 fout << fMateIDs[ic + ir * fNoVoxelsX /
99 if (ic != fNoVoxelsX / fCompress - 1) fo
100 }
101 fout << G4endl;
102 }
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oo
106 void DicomFileCT::DumpDensitiesToTextFile(std:
107 {
108 G4int fCompress = theFileMgr->GetCompression
109 if (DicomFileMgr::verbose >= warningVerb)
110 G4cout << fLocation << " DumpDensitiesToTe
111 << G4endl;
112
113 G4int copyNo = 0;
114 for (int ir = 0; ir < fNoVoxelsY / fCompress
115 for (int ic = 0; ic < fNoVoxelsX / fCompre
116 fout << fDensities[ic + ir * fNoVoxelsX
117 if (ic != fNoVoxelsX / fCompress - 1) fo
118 if (copyNo % 8 == 7) fout << G4endl;
119 copyNo++;
120 }
121 if (copyNo % 8 != 0) fout << G4endl;
122 }
123 }
124
125 //....oooOO0OOooo........oooOO0OOooo........oo
126 void DicomFileCT::BuildStructureIDs()
127 {
128 G4int fCompress = theFileMgr->GetCompression
129 std::vector<DicomFileStructure*> dfs = theFi
130 if (dfs.size() == 0) return;
131
132 G4int NMAXROI = DicomFileMgr::GetInstance()-
133 G4int NMAXROI_real = std::log(INT_MAX) / std
134
135 // fStructure = new G4int(fNoVoxelsX*fNoVox
136 for (int ir = 0; ir < fNoVoxelsY / fCompress
137 for (int ic = 0; ic < fNoVoxelsX / fCompre
138 // fStructure[ic+ir*fNoVoxelsX] = 0
139 fStructure.push_back(0);
140 }
141 }
142
143 std::set<double> distInters;
144
145 // std::fill_n(fStructure,fNoVoxelsX*fNoVox
146 //
147 for (size_t ii = 0; ii < dfs.size(); ii++) {
148 std::vector<DicomROI*> rois = dfs[ii]->Get
149 for (size_t jj = 0; jj < rois.size(); jj++
150 if (DicomFileMgr::verbose >= debugVerb)
151 std::cout << " BuildStructureIDs check
152 << rois[jj]->GetName() << G4
153 std::vector<DicomROIContour*> contours =
154 // G4int roi = jj+1;
155 G4int roiID = rois[jj]->GetNumber();
156 for (size_t kk = 0; kk < contours.size()
157 // check contour corresponds to this C
158 DicomROIContour* roic = contours[kk];
159 // if( DicomVerb(-debugVerb) ) G4cout
160 // << " " << fLocation
161 if (DicomFileMgr::verbose >= debugVerb
162 G4cout << " " << roiID << " " << kk
163 << fMinZ << " " << fMaxZ << G
164 // Check Contour correspond to slice
165
166 if (roic->GetZ() > fMaxZ || roic->GetZ
167 if (DicomFileMgr::verbose >= debugVerb
168 G4cout << " INTERS CONTOUR WITH Z SL
169 << fMaxZ << G4endl;
170 if (roic->GetGeomType() == "CLOSED_PLA
171 // Get min and max X and Y, and corr
172 std::vector<G4ThreeVector> points =
173 if (DicomFileMgr::verbose >= debugVe
174 G4cout << jj << " " << kk << " NPO
175 std::vector<G4ThreeVector> dirs = ro
176 double minXc = DBL_MAX;
177 double maxXc = -DBL_MAX;
178 double minYc = DBL_MAX;
179 double maxYc = -DBL_MAX;
180 for (size_t ll = 0; ll < points.size
181 minXc = std::min(minXc, points[ll]
182 maxXc = std::max(maxXc, points[ll]
183 minYc = std::min(minYc, points[ll]
184 maxYc = std::max(maxYc, points[ll]
185 }
186 if (minXc < fMinX || maxXc > fMaxX |
187 G4cerr << " minXc " << minXc << "
188 << " minYc " << minYc << "
189 << G4endl;
190 G4Exception("DicomFileCT::BuildStr
191 "Contour limits exceed
192 }
193 int idMinX = std::max(0, int((minXc
194 int idMaxX =
195 std::min(fNoVoxelsX / fCompress -
196 int idMinY = std::max(0, int((minYc
197 int idMaxY =
198 std::min(fNoVoxelsY / fCompress -
199 if (DicomFileMgr::verbose >= debugVe
200 G4cout << " minXc " << minXc << "
201 << " minYc " << minYc << "
202 << G4endl;
203 if (DicomFileMgr::verbose >= debugVe
204 G4cout << " idMinX " << idMinX <<
205 << " idMaxY " << idMaxY <<
206 // for each voxel: build 4 lines fro
207 // and check how many contour segme
208 for (int ix = idMinX; ix <= idMaxX;
209 for (int iy = idMinY; iy <= idMaxY
210 G4bool bOK = false;
211 G4int bOKs;
212 if (DicomFileMgr::verbose >= deb
213 G4cout << "@@@@@ TRYING POINT
214 << "," << fMinY + fVoxe
215 // four corners
216 for (int icx = 0; icx <= 1; icx+
217 for (int icy = 0; icy <= 1; ic
218 bOKs = 0;
219 if (bOK) continue;
220 double p0x = fMinX + fVoxelD
221 double p0y = fMinY + fVoxelD
222 double v0x = 1.;
223 if (icx == 1) v0x = -1.;
224 double v0y = 0.99 * fVoxelDi
225 if (DicomFileMgr::verbose >=
226 G4cout << ix << " + " << i
227 << "," << p0y << ")
228 << " DIR= (" << v0x
229 int NTRIES = theFileMgr->Get
230 for (int ip = 0; ip < NTRIES
231 distInters.clear();
232 v0y -= ip * 0.1 * v0y;
233 G4double halfDiagonal = 0.
234 if (DicomFileMgr::verbose
235 G4cout << ip << " TRYING
236 << " DIR= (" << v
237 for (size_t ll = 0; ll < p
238 double d0x = points[ll].
239 double d0y = points[ll].
240 double w0x = dirs[ll].x(
241 double w0y = dirs[ll].y(
242 double fac1 = w0x * v0y
243 if (fac1 == 0) { // par
244 continue;
245 }
246 double fac2 = d0x * v0y
247 double fac3 = d0y * w0x
248 double lambdaq = -fac2 /
249 if (lambdaq < 0. || lamb
250 // intersection further
251 double lambdap = fac3 /
252 if (lambdap > 0.) {
253 distInters.insert(lamb
254 if (DicomFileMgr::verb
255 G4cout << " !! GOOD
256 << d0x + p0x
257 << ") = ("
258 << ") "
259 << " N " << d
260 }
261 if (DicomFileMgr::verbos
262 G4cout << " INTERS LAM
263
264 if (DicomFileMgr::verbos
265 G4cout << " INTERS POI
266 << d0y + p0y +
267 << "," << p0y +
268 }
269 if (distInters.size() % 2
270 if (*(distInters.begin()
271 //
272 bOKs += 1;
273 if (DicomFileMgr::verb
274 G4cout << " OK= " <<
275 << *(distInte
276 }
277 }
278 }
279 if (bOKs == NTRIES) {
280 bOK = true;
281 if (DicomFileMgr::verbose
282 G4cout << "@@@@@ POINT O
283 << "," << fMinY +
284 }
285 }
286 } // loop to four corners
287 if (bOK) {
288 // extract previous ROI value
289 int roival = fStructure[ix + i
290 // roival = 2 +
291 if (roival != 0 && roival != i
292 std::set<G4int> roisVoxel;
293 int nRois = std::log10(roiva
294 if (nRois > NMAXROI_real) {
295 G4Exception("DicomFileCT:B
296 G4String("Too
297 " + std::to_string(nRois) + " > " + std::to_st
298 + " T
299 GING -NStructureNMaxROI argument to a lower va
300 .c_str());
301 }
302 for (int inr = 0; inr < nRoi
303 roisVoxel.insert(int(roiva
304 }
305 roisVoxel.insert(roiID);
306 roival = 0;
307 size_t inr = 0;
308 for (std::set<G4int>::const_
309 ite != roisVoxel.end();
310 {
311 roival += (*ite) * std::po
312 }
313 fStructure[ix + iy * fNoVoxe
314 if (DicomFileMgr::verbose >=
315 G4cout << " WITH PREVIOUS
316 }
317 }
318 else {
319 fStructure[ix + iy * fNoVoxe
320 }
321 }
322 }
323 }
324 }
325 }
326 }
327 }
328
329 if (DicomFileMgr::verbose >= 1) {
330 //@@@@ PRINT structures
331 //@@@ PRINT points of structures in this Z
332 if (DicomFileMgr::verbose >= 0) G4cout <<
333 for (size_t ii = 0; ii < dfs.size(); ii++)
334 std::vector<DicomROI*> rois = dfs[ii]->G
335 for (size_t jj = 0; jj < rois.size(); jj
336 int roi = rois[jj]->GetNumber();
337 std::vector<DicomROIContour*> contours
338 for (size_t kk = 0; kk < contours.size
339 DicomROIContour* roic = contours[kk]
340 // check contour corresponds to this
341 if (roic->GetZ() > fMaxZ || roic->Ge
342 if (roic->GetGeomType() == "CLOSED_P
343 if (DicomFileMgr::verbose >= testV
344 G4cout << " PRINTING CONTOUR? "
345 << roic->GetZ() << " SLIC
346 if (roic->GetZ() > fMaxZ || roic->
347 std::vector<G4ThreeVector> points
348 std::vector<G4ThreeVector> dirs =
349 if (DicomFileMgr::verbose >= 0)
350 std::cout << " STRUCTURE Z " <<
351 for (size_t ll = 0; ll < points.si
352 if (DicomFileMgr::verbose >= 0)
353 G4cout << roi << " : " << ll <
354 << points[ll].y() << ")
355 << " (" << dirs[ll].x()
356 }
357 }
358 }
359 }
360 }
361 //@@@ PRINT points in slice inside structu
362 for (int ir = 0; ir < fNoVoxelsY / fCompre
363 for (int ic = 0; ic < fNoVoxelsX / fComp
364 if (fStructure[ic + ir * fNoVoxelsX /
365 if (DicomFileMgr::verbose >= 0)
366 G4cout << ic + ir * fNoVoxelsX / f
367 << " STRUCTURE VOXEL (" <<
368 << fMinY + fVoxelDimY * fCo
369 << ") = " << fStructure[ic
370 }
371 }
372 }
373 }
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oo
377 void DicomFileCT::DumpStructureIDsToTextFile(s
378 {
379 G4int fCompress = theFileMgr->GetCompression
380 if (DicomFileMgr::verbose >= 0)
381 G4cout << fLocation << " DumpStructureIDsT
382 << G4endl;
383 std::vector<DicomFileStructure*> dfs = theFi
384 if (dfs.size() == 0) return;
385
386 for (int ir = 0; ir < fNoVoxelsY / fCompress
387 for (int ic = 0; ic < fNoVoxelsX / fCompre
388 fout << fStructure[ic + ir * fNoVoxelsX
389 if (ic != fNoVoxelsX / fCompress - 1) fo
390 }
391 fout << G4endl;
392 }
393 }
394