| Geant4 Cross Reference |
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 "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........oooOO0OOooo........oooOO0OOooo......
42 DicomFileCT::DicomFileCT() {}
43
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 DicomFileCT::DicomFileCT(DcmDataset* dset) : DicomVFileImage(dset) {}
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 void DicomFileCT::BuildMaterials()
49 {
50 G4int fCompress = theFileMgr->GetCompression();
51 if (fNoVoxelsX % fCompress != 0 || fNoVoxelsY % fCompress != 0) {
52 G4Exception("DicompFileMgr.:BuildMaterials", "DFC004", FatalException,
53 ("Compression factor = " + std::to_string(fCompress)
54 + " has to be a divisor of Number of voxels X = " + std::to_string(fNoVoxelsX)
55 + " and Y " + std::to_string(fNoVoxelsY))
56 .c_str());
57 }
58
59 // if( DicomVerb(debugVerb) ) G4cout << " BuildMaterials " << fFileName << G4endl;
60 double meanHV = 0.;
61 for (int ir = 0; ir < fNoVoxelsY; ir += fCompress) {
62 for (int ic = 0; ic < fNoVoxelsX; ic += fCompress) {
63 meanHV = 0.;
64 int isumrMax = std::min(ir + fCompress, fNoVoxelsY);
65 int isumcMax = std::min(ic + fCompress, fNoVoxelsX);
66 for (int isumr = ir; isumr < isumrMax; isumr++) {
67 for (int isumc = ic; isumc < isumcMax; isumc++) {
68 meanHV += fHounsfieldV[isumc + isumr * fNoVoxelsX];
69 // G4cout << isumr << " " << isumc << " GET mean " << meanHV << G4endl;
70 }
71 }
72 meanHV /= (isumrMax - ir) * (isumcMax - ic);
73 G4double meanDens = theFileMgr->Hounsfield2density(meanHV);
74 // G4cout << ir << " " << ic << " FINAL mean " << meanDens << G4endl;
75
76 fDensities.push_back(meanDens);
77 size_t mateID;
78 if (theFileMgr->IsMaterialsDensity()) {
79 mateID = theFileMgr->GetMaterialIndexByDensity(meanDens);
80 }
81 else {
82 mateID = theFileMgr->GetMaterialIndex(meanHV);
83 }
84 fMateIDs.push_back(mateID);
85 }
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 void DicomFileCT::DumpMateIDsToTextFile(std::ofstream& fout)
91 {
92 G4int fCompress = theFileMgr->GetCompression();
93 if (DicomFileMgr::verbose >= warningVerb)
94 G4cout << fLocation << " DumpMateIDsToTextFile " << fFileName << " " << fMateIDs.size()
95 << G4endl;
96 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
97 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
98 fout << fMateIDs[ic + ir * fNoVoxelsX / fCompress];
99 if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
100 }
101 fout << G4endl;
102 }
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 void DicomFileCT::DumpDensitiesToTextFile(std::ofstream& fout)
107 {
108 G4int fCompress = theFileMgr->GetCompression();
109 if (DicomFileMgr::verbose >= warningVerb)
110 G4cout << fLocation << " DumpDensitiesToTextFile " << fFileName << " " << fDensities.size()
111 << G4endl;
112
113 G4int copyNo = 0;
114 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
115 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
116 fout << fDensities[ic + ir * fNoVoxelsX / fCompress];
117 if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
118 if (copyNo % 8 == 7) fout << G4endl;
119 copyNo++;
120 }
121 if (copyNo % 8 != 0) fout << G4endl;
122 }
123 }
124
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126 void DicomFileCT::BuildStructureIDs()
127 {
128 G4int fCompress = theFileMgr->GetCompression();
129 std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
130 if (dfs.size() == 0) return;
131
132 G4int NMAXROI = DicomFileMgr::GetInstance()->GetStructureNMaxROI();
133 G4int NMAXROI_real = std::log(INT_MAX) / std::log(NMAXROI);
134
135 // fStructure = new G4int(fNoVoxelsX*fNoVoxelsY);
136 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
137 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
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*fNoVoxelsY,0);
146 //
147 for (size_t ii = 0; ii < dfs.size(); ii++) {
148 std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
149 for (size_t jj = 0; jj < rois.size(); jj++) {
150 if (DicomFileMgr::verbose >= debugVerb)
151 std::cout << " BuildStructureIDs checking ROI " << rois[jj]->GetNumber() << " "
152 << rois[jj]->GetName() << G4endl;
153 std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
154 // G4int roi = jj+1;
155 G4int roiID = rois[jj]->GetNumber();
156 for (size_t kk = 0; kk < contours.size(); kk++) {
157 // check contour corresponds to this CT slice
158 DicomROIContour* roic = contours[kk];
159 // if( DicomVerb(-debugVerb) ) G4cout << jj << " " << kk << " INTERS CONTOUR " << roic
160 // << " " << fLocation << G4endl;
161 if (DicomFileMgr::verbose >= debugVerb)
162 G4cout << " " << roiID << " " << kk << " INTERS CONTOUR " << roic->GetZ() << " SLICE Z "
163 << fMinZ << " " << fMaxZ << G4endl;
164 // Check Contour correspond to slice
165
166 if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
167 if (DicomFileMgr::verbose >= debugVerb)
168 G4cout << " INTERS CONTOUR WITH Z SLIZE " << fMinZ << " < " << roic->GetZ() << " < "
169 << fMaxZ << G4endl;
170 if (roic->GetGeomType() == "CLOSED_PLANAR") {
171 // Get min and max X and Y, and corresponding slice indexes
172 std::vector<G4ThreeVector> points = roic->GetPoints();
173 if (DicomFileMgr::verbose >= debugVerb)
174 G4cout << jj << " " << kk << " NPOINTS " << points.size() << G4endl;
175 std::vector<G4ThreeVector> dirs = roic->GetDirections();
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(); ll++) {
181 minXc = std::min(minXc, points[ll].x());
182 maxXc = std::max(maxXc, points[ll].x());
183 minYc = std::min(minYc, points[ll].y());
184 maxYc = std::max(maxYc, points[ll].y());
185 }
186 if (minXc < fMinX || maxXc > fMaxX || minYc < fMinY || maxYc > fMaxY) {
187 G4cerr << " minXc " << minXc << " < " << fMinX << " maxXc " << maxXc << " > " << fMaxX
188 << " minYc " << minYc << " < " << fMinY << " maxYc " << maxYc << " > " << fMaxY
189 << G4endl;
190 G4Exception("DicomFileCT::BuildStructureIDs", "DFCT001", JustWarning,
191 "Contour limits exceed Z slice extent");
192 }
193 int idMinX = std::max(0, int((minXc - fMinX) / fVoxelDimX / fCompress));
194 int idMaxX =
195 std::min(fNoVoxelsX / fCompress - 1, int((maxXc - fMinX) / fVoxelDimX / fCompress + 1));
196 int idMinY = std::max(0, int((minYc - fMinY) / fVoxelDimY / fCompress));
197 int idMaxY =
198 std::min(fNoVoxelsY / fCompress - 1, int((maxYc - fMinY) / fVoxelDimY / fCompress + 1));
199 if (DicomFileMgr::verbose >= debugVerb)
200 G4cout << " minXc " << minXc << " < " << fMinX << " maxXc " << maxXc << " > " << fMaxX
201 << " minYc " << minYc << " < " << fMinY << " maxYc " << maxYc << " > " << fMaxY
202 << G4endl;
203 if (DicomFileMgr::verbose >= debugVerb)
204 G4cout << " idMinX " << idMinX << " idMaxX " << idMaxX << " idMinY " << idMinY
205 << " idMaxY " << idMaxY << G4endl;
206 // for each voxel: build 4 lines from the corner towards the center
207 // and check how many contour segments it crosses, and the minimum distance to a segment
208 for (int ix = idMinX; ix <= idMaxX; ix++) {
209 for (int iy = idMinY; iy <= idMaxY; iy++) {
210 G4bool bOK = false;
211 G4int bOKs;
212 if (DicomFileMgr::verbose >= debugVerb)
213 G4cout << "@@@@@ TRYING POINT (" << fMinX + fVoxelDimX * fCompress * (ix + 0.5)
214 << "," << fMinY + fVoxelDimY * fCompress * (iy + 0.5) << ")" << G4endl;
215 // four corners
216 for (int icx = 0; icx <= 1; icx++) {
217 for (int icy = 0; icy <= 1; icy++) {
218 bOKs = 0;
219 if (bOK) continue;
220 double p0x = fMinX + fVoxelDimX * fCompress * (ix + icx);
221 double p0y = fMinY + fVoxelDimY * fCompress * (iy + icy);
222 double v0x = 1.;
223 if (icx == 1) v0x = -1.;
224 double v0y = 0.99 * fVoxelDimY / fVoxelDimX * std::pow(-1., icy);
225 if (DicomFileMgr::verbose >= testVerb)
226 G4cout << ix << " + " << icx << " " << iy << " + " << icy << " CORNER (" << p0x
227 << "," << p0y << ") "
228 << " DIR= (" << v0x << "," << v0y << ") " << G4endl;
229 int NTRIES = theFileMgr->GetStructureNCheck();
230 for (int ip = 0; ip < NTRIES; ip++) {
231 distInters.clear();
232 v0y -= ip * 0.1 * v0y;
233 G4double halfDiagonal = 0.5 * fVoxelDimX * fCompress;
234 if (DicomFileMgr::verbose >= testVerb)
235 G4cout << ip << " TRYING WITH DIRECTION ("
236 << " DIR= (" << v0x << "," << v0y << ") " << bOKs << G4endl;
237 for (size_t ll = 0; ll < points.size(); ll++) {
238 double d0x = points[ll].x() - p0x;
239 double d0y = points[ll].y() - p0y;
240 double w0x = dirs[ll].x();
241 double w0y = dirs[ll].y();
242 double fac1 = w0x * v0y - w0y * v0x;
243 if (fac1 == 0) { // parallel lines
244 continue;
245 }
246 double fac2 = d0x * v0y - d0y * v0x;
247 double fac3 = d0y * w0x - d0x * w0y;
248 double lambdaq = -fac2 / fac1;
249 if (lambdaq < 0. || lambdaq >= 1.) continue;
250 // intersection further than segment length
251 double lambdap = fac3 / fac1;
252 if (lambdap > 0.) {
253 distInters.insert(lambdap);
254 if (DicomFileMgr::verbose >= testVerb)
255 G4cout << " !! GOOD INTERS " << lambdaq << " ("
256 << d0x + p0x + lambdaq * w0x << "," << d0y + p0y + lambdaq * w0y
257 << ") = (" << p0x + lambdap * v0x << "," << p0y + lambdap * v0y
258 << ") "
259 << " N " << distInters.size() << G4endl;
260 }
261 if (DicomFileMgr::verbose >= testVerb)
262 G4cout << " INTERS LAMBDAQ " << lambdaq << " P " << lambdap << G4endl;
263
264 if (DicomFileMgr::verbose >= debugVerb)
265 G4cout << " INTERS POINT (" << d0x + p0x + lambdaq * w0x << ","
266 << d0y + p0y + lambdaq * w0y << ") = (" << p0x + lambdap * v0x
267 << "," << p0y + lambdap * v0y << ") " << G4endl;
268 }
269 if (distInters.size() % 2 == 1) {
270 if (*(distInters.begin()) > halfDiagonal) {
271 // bOK = true;
272 bOKs += 1;
273 if (DicomFileMgr::verbose >= debugVerb)
274 G4cout << " OK= " << bOKs << " N INTERS " << distInters.size() << " "
275 << *(distInters.begin()) << " > " << halfDiagonal << G4endl;
276 }
277 }
278 }
279 if (bOKs == NTRIES) {
280 bOK = true;
281 if (DicomFileMgr::verbose >= debugVerb)
282 G4cout << "@@@@@ POINT OK (" << fMinX + fVoxelDimX * fCompress * (ix + 0.5)
283 << "," << fMinY + fVoxelDimY * fCompress * (iy + 0.5) << ")" << G4endl;
284 }
285 }
286 } // loop to four corners
287 if (bOK) {
288 // extract previous ROI value
289 int roival = fStructure[ix + iy * fNoVoxelsX / fCompress];
290 // roival = 2 + NMAXROI*3 + NMAXROI*NMAXROI*15;
291 if (roival != 0 && roival != int(roiID)) {
292 std::set<G4int> roisVoxel;
293 int nRois = std::log10(roival) / std::log10(NMAXROI) + 1;
294 if (nRois > NMAXROI_real) { // another one cannot be added
295 G4Exception("DicomFileCT:BuildStructureIDs", "DFC0004", FatalException,
296 G4String("Too many ROIs associated to a voxel: \
297 " + std::to_string(nRois) + " > " + std::to_string(NMAXROI_real)
298 + " TRY CHAN\
299 GING -NStructureNMaxROI argument to a lower value")
300 .c_str());
301 }
302 for (int inr = 0; inr < nRois; inr++) {
303 roisVoxel.insert(int(roival / std::pow(NMAXROI, inr)) % NMAXROI);
304 }
305 roisVoxel.insert(roiID);
306 roival = 0;
307 size_t inr = 0;
308 for (std::set<G4int>::const_iterator ite = roisVoxel.begin();
309 ite != roisVoxel.end(); ite++, inr++)
310 {
311 roival += (*ite) * std::pow(NMAXROI, inr);
312 }
313 fStructure[ix + iy * fNoVoxelsX / fCompress] = roival;
314 if (DicomFileMgr::verbose >= testVerb) {
315 G4cout << " WITH PREVIOUS ROI IN VOXEL " << roival << G4endl;
316 }
317 }
318 else {
319 fStructure[ix + iy * fNoVoxelsX / fCompress] = roiID;
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 slice
332 if (DicomFileMgr::verbose >= 0) G4cout << " STRUCTURES FOR SLICE " << fLocation << G4endl;
333 for (size_t ii = 0; ii < dfs.size(); ii++) {
334 std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
335 for (size_t jj = 0; jj < rois.size(); jj++) {
336 int roi = rois[jj]->GetNumber();
337 std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
338 for (size_t kk = 0; kk < contours.size(); kk++) {
339 DicomROIContour* roic = contours[kk];
340 // check contour corresponds to this CT slice
341 if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
342 if (roic->GetGeomType() == "CLOSED_PLANAR") {
343 if (DicomFileMgr::verbose >= testVerb)
344 G4cout << " PRINTING CONTOUR? " << roi << " " << kk << " INTERS CONTOUR "
345 << roic->GetZ() << " SLICE Z " << fMinZ << " " << fMaxZ << G4endl;
346 if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
347 std::vector<G4ThreeVector> points = roic->GetPoints();
348 std::vector<G4ThreeVector> dirs = roic->GetDirections();
349 if (DicomFileMgr::verbose >= 0)
350 std::cout << " STRUCTURE Z " << roic->GetZ() << std::endl;
351 for (size_t ll = 0; ll < points.size(); ll++) {
352 if (DicomFileMgr::verbose >= 0)
353 G4cout << roi << " : " << ll << " STRUCTURE POINT (" << points[ll].x() << ","
354 << points[ll].y() << ") "
355 << " (" << dirs[ll].x() << "," << dirs[ll].y() << ") = " << roi << G4endl;
356 }
357 }
358 }
359 }
360 }
361 //@@@ PRINT points in slice inside structure
362 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
363 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
364 if (fStructure[ic + ir * fNoVoxelsX / fCompress] != 0) {
365 if (DicomFileMgr::verbose >= 0)
366 G4cout << ic + ir * fNoVoxelsX / fCompress << " = " << ic << " " << ir
367 << " STRUCTURE VOXEL (" << fMinX + fVoxelDimX * fCompress * (ic + 0.5) << ","
368 << fMinY + fVoxelDimY * fCompress * (ir + 0.5)
369 << ") = " << fStructure[ic + ir * fNoVoxelsX / fCompress] << G4endl;
370 }
371 }
372 }
373 }
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 void DicomFileCT::DumpStructureIDsToTextFile(std::ofstream& fout)
378 {
379 G4int fCompress = theFileMgr->GetCompression();
380 if (DicomFileMgr::verbose >= 0)
381 G4cout << fLocation << " DumpStructureIDsToTextFile " << fFileName << " " << fStructure.size()
382 << G4endl;
383 std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
384 if (dfs.size() == 0) return;
385
386 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
387 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
388 fout << fStructure[ic + ir * fNoVoxelsX / fCompress];
389 if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
390 }
391 fout << G4endl;
392 }
393 }
394